{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Many Particle Model (MPM) "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Many Paticle Model (MPM) of a lithium-ion battery is an extension of the Single Particle Model to account for a continuous distribution of active particle sizes in each electrode $\\text{k}=\\text{n},\\text{p}$. Therefore, many of the same model assumptions hold, e.g., the transport in the electrolyte is instantaneous and hence the through-cell variation (in $x$) is neglected. The full set of assumptions and description of the particle size geometry is given in [[4]](#References). Note that the MPM in [[4]](#References) is for a half cell and the version implemented in PyBaMM is for a full cell and uses the notation and scaling given in [[5]](#References).\n",
    "\n",
    "\n",
    "## Particle size geometry\n",
    "\n",
    "In each electrode $\\text{k}=\\text{n},\\text{p}$, there are spherical particles of each radius $R_\\text{k}$ in the range $R_\\text{k,min}<R_\\text{k}<R_\\text{k,max}$, with the fraction of all particles of a given radius $R_\\text{k}$ given by the particle-size distribution (base on number)\n",
    "$f_\\text{k,num}(R_\\text{k})$. However, it is more convenient to deal with the fraction\n",
    "of _surface area_ contributed by particles of radius $R_\\text{k}$, which we denote $f_{\\text{k},a}(R_\\text{k})$ and refer to as the _area-weighted_ particle-size distribution. The area and number-based distributions are related via\n",
    "$$\n",
    "f_{\\text{k},a}(R_\\text{k}) = \\frac{4\\pi R_\\text{k}^2 f_\\text{k,num}(R_\\text{k})}{\\int_{R_\\text{k,min}}^{R_\\text{k,max}} 4\\pi R_\\text{k}^2 f_\\text{k,num}(R_\\text{k})\\,\\text{d}R_\\text{k}}\n",
    "$$\n",
    "The total amount of surface area (per unit volume) $a_\\text{k,tot}$ accounting for all particles is expressed in terms of the active material volume fraction $\\epsilon_{s,\\text{k}}$, similar to the other models in PyBaMM (SPM, DFN):\n",
    "$$\n",
    "\\epsilon_{s,\\text{k}}= \\int \\frac{1}{3} R_\\text{k} \\underbrace{a_\\text{k,tot}f_{\\text{k},a}(R_\\text{k})}_{\\text{area }a_\\text{k}(R_\\text{k})\\text{ of particles size }R_\\text{k}}\\,\\text{d}R_\\text{k}\n",
    "$$\n",
    "Rearranging and defining $\\bar{R}_{\\text{k},a}=\\int  R_\\text{k} f_{\\text{k},a}(R_\\text{k})\\,\\text{d}R_\\text{k}$ as the mean of the area-weighted distribution, we find\n",
    "$$\n",
    "a_\\text{k,tot}=\\frac{3\\epsilon_{s,\\text{k}}}{\\int  R_\\text{k} f_{\\text{k},a}(R_\\text{k})\\,\\text{d}R_\\text{k}} = \\frac{3\\epsilon_{s,\\text{k}}}{\\bar{R}_{\\text{k},a}}.\n",
    "$$\n",
    "Then $a_\\text{k,tot}$ is the aggregate surface area of the particle population and analogous to the variables `\"X-averaged negative electrode surface area to volume ratio [m-1]\"`, etc. in the SPM, SPMe, and DFN models, and can be calculated in a similar way as shown above using the _area-weighted mean radius_ $\\bar{R}_{\\text{k},a}$ (other mean radii do not have this property). See [[4]](#References) for more details on the different types of distribution and mean radii. \n",
    "\n",
    "Another common way to express the size distribution is via particle volume. The fraction of volume contributed by the particles of radius $R_\\text{k}$, denoted the _volume-weighted_ particle-size distribution is related to the number and area ones via\n",
    "$$\n",
    "f_{\\text{k},v}(R_\\text{k}) = \\frac{\\frac{1}{3} R_\\text{k} f_{\\text{k},a}(R_\\text{k})}{\\int_{R_\\text{k,min}}^{R_\\text{k,max}} \\frac{1}{3} R_\\text{k} f_{\\text{k},a}(R_\\text{k})\\,\\text{d}R_\\text{k}} =\\frac{\\frac{4}{3}\\pi R_\\text{k}^3 f_\\text{k,num}(R_\\text{k})}{\\int_{R_\\text{k,min}}^{R_\\text{k,max}} \\frac{4}{3}\\pi R_\\text{k}^3 f_\\text{k,num}(R_\\text{k})\\,\\text{d}R_\\text{k}}\n",
    "$$\n",
    "\n",
    "It is sufficient to specify $f_{\\text{k},a}(R_\\text{k})$, which is present requirement in the MPM."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model equations\n",
    "\n",
    "In each electrode, only one representative particle of each size $R_\\text{k}$ needs to be modelled. The concentration of lithium in the solid particles is denoted $c_{\\text{s,k}}(t,R_\\text{k}, r_\\text{k})$, which varies with time, particle radius $R_\\text{k}$, and the radial coordinate $r_{\\text{k}} \\in[0,R_{\\text{k}}]$ within the spherical particle. The potential is uniform across all particles in the electrode, $\\phi_{\\text{s,k}}(t)$. \n",
    "\n",
    "The equations for molar conservation of lithium ($c_{\\text{s,k}}$) are then: \n",
    "$$\n",
    "\\frac{\\partial c_{\\text{s,k}}}{\\partial t} = -\\frac{1}{r_{\\text{k}}^2} \\frac{\\partial}{\\partial r_{\\text{k}}} \\left(r_{\\text{k}}^2 N_{\\text{s,k}}\\right), \\\\\n",
    "N_{\\text{s,k}} = -D_{\\text{s,k}}(c_{\\text{s,k}}) \\frac{\\partial c_{\\text{s,k}}}{\\partial r_{\\text{k}}}, \\quad \\text{k} \\in \\text{n, p}, \n",
    "$$\n",
    "$$\n",
    "N_{\\text{s,k}}\\big|_{r_{\\text{k}}=0} = 0,  \\ \\  N_{\\text{s,k}}\\big|_{r_{\\text{k}}=R_{\\text{k}}} = \n",
    "\\frac{j_{\\text{k}}}{F} \\quad \\text{k} \\in \\text{n, p}, \\quad\\\\\n",
    "c_{\\text{s,k}}(0,R_\\text{k},r_{\\text{k}}) = c_{\\text{s,k,0}}, \\quad \\text{k} \\in \\text{n, p},$$\n",
    "where $D_{\\text{s,k}}$ is the diffusion coefficient in the solid, $N_{\\text{s,k}}$ denotes the flux of lithium ions in the solid particle, $F$ is Faraday's constant. The interfacial current density is given by $j_\\text{k}$, which also varies with particle size.\n",
    "\n",
    "### Algebraic equations for the potentials\n",
    "The potentials $\\phi_{\\text{s,k}}(t)$ are determined via the integral constraint that the total current flowing across the electrode interface must equal (up to a minus sign) the through-cell current density $i$. Writing this in terms of the potential differences $\\Delta \\phi_{\\text{s,k}} = \\phi_{\\text{s,k}} - \\phi_{\\text{e}}$,\n",
    "$$\n",
    "L_\\text{k}a_\\text{k,tot}\\int_{R_\\text{k,min}}^{R_\\text{k,max}} f_{\\text{k},a}(R_\\text{k})j_\\text{k}\\,\\text{d}R_\\text{k} = \\begin{cases}\n",
    "i,\\quad \\text{k}=\\text{n}\\\\\n",
    "-i,\\quad \\text{k}=\\text{p}\n",
    "\\end{cases}\n",
    "$$\n",
    "with Butler-Volmer kinetics\n",
    "$$\n",
    "j_\\text{k}=j_{\\text{0,k}} \\sinh\\left[\\frac{F}{2R_g T}(\\Delta \\phi_{\\text{s,k}}-U_{\\text{k}}(c_{\\text{s},\\text{k}}))\\right], \\ \\ j_{\\text{0,k}} =  m_{\\text{k}}(c_{\\text{e}}c_{\\text{s,k}})^{1/2}(c_\\text{k,max}-c_{\\text{s,k}})^{1/2}.\n",
    "$$\n",
    "This gives an integral (or algebraic once discretized) equation for $\\Delta \\phi_{\\text{s,k}}$ which is coupled to the concentration equations above.\n",
    "The voltage is then obtained from\n",
    "$$\n",
    "V = \\Delta \\phi_{\\text{s,p}} - \\Delta \\phi_{\\text{s,n}}\n",
    "$$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example solving MPM\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Note: you may need to restart the kernel to use updated packages.\n"
     ]
    }
   ],
   "source": [
    "%pip install \"pybamm[plot,cite]\" -q    # install PyBaMM if it is not installed\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import pybamm"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create an instance of the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = pybamm.lithium_ion.MPM()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's inspect some variables (e.g. the lithium concentration and interfacial current densities) that depend on particle size $R_\\text{k}$. The variables of interest are `X-averaged` versions as there is no dependence on $x$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X-averaged negative particle concentration\n",
      "X-averaged negative particle concentration [mol.m-3]\n",
      "X-averaged negative particle concentration distribution\n",
      "X-averaged negative particle concentration distribution [mol.m-3]\n"
     ]
    }
   ],
   "source": [
    "model.variables.search(\"X-averaged negative particle concentration\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The concentration that is being solved for in the MPM, and which varies with particle size is the one ending in `\"distribution\"`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'primary': ['negative particle'],\n",
       " 'secondary': ['negative particle size'],\n",
       " 'tertiary': ['current collector']}"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "c_n_R_dependent = model.variables[\n",
    "    \"X-averaged negative particle concentration distribution [mol.m-3]\"\n",
    "]\n",
    "c_n_R_dependent.domains"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that the secondary domain is `'negative particle size'`, which is treated as another (microscale) domain in PyBaMM.\n",
    "\n",
    "The variable without the `\"distribution\"` has been \"size averaged\" and can be compared to the variable with the same name from the other lithium-ion models in PyBaMM with only a single particle size. The concentration within the particles is a volume-based quantity and is thus averaged by volume (to preserve the total amount of lithium):\n",
    "$$\n",
    "\\left<c_{\\text{s,k}}\\right>_v = \\int_{R_\\text{k,min}}^{R_\\text{k,max}} f_{\\text{k},v}(R_\\text{k})c_{\\text{s,k}}(t,R_\\text{k}, r_\\text{k})\\,\\text{d}R_\\text{k}\n",
    "$$\n",
    "\n",
    "In particular, if the variance of the particle-size distribution $f_{\\text{k},a}$ is shrunk to zero and all particles become concentrated at its mean radius $\\bar{R}_{\\text{k},a}$, the variable `\"X-averaged negative particle concentration [mol.m-3]\"` will coincide with the same variable from an SPM with particle radius $R_\\text{k}=\\bar{R}_{\\text{k},a}$. However, `\"X-averaged negative particle concentration distribution [mol.m-3]\"` will remain \"particle-size dependent\"."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The convention of adding `\"distribution\"` to the end of a variable name to indicate particle-size dependence has been used for other variables, such as the interfacial current density:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X-averaged negative electrode interfacial current density\n",
      "X-averaged negative electrode interfacial current density [A.m-2]\n",
      "X-averaged negative electrode interfacial current density distribution\n",
      "X-averaged negative electrode interfacial current density distribution [A.m-2]\n",
      "X-averaged negative electrode interfacial current density per volume [A.m-3]\n"
     ]
    }
   ],
   "source": [
    "model.variables.search(\"X-averaged negative electrode interfacial current density\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As the interfacial current density is a flux per unit area on the particle surface, the \"size averaging\" is done by area (to preserve the total flux of lithium):\n",
    "$$\n",
    "\\left<j_{\\text{k}}\\right>_a = \\int_{R_\\text{k,min}}^{R_\\text{k,max}} f_{\\text{k},a}(R_\\text{k})j_{\\text{k}}(t,R_\\text{k})\\,\\text{d}R_\\text{k}\n",
    "$$\n",
    "The averaging is merely done to allow comparison to variables from other models with only a single size, and are not necessarily used within the MPM itself, or are physically meaningful.\n",
    "\n",
    "Note: not all variables have a \"distribution\" version, such as the potentials or temperature variables, as they do not vary with particle size in the MPM as implemented here."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mesh points\n",
    " By default, the size domain is discretized into 30 grid points on a uniform 1D mesh."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "negative electrode is of type Uniform1DSubMesh\n",
      "separator is of type Uniform1DSubMesh\n",
      "positive electrode is of type Uniform1DSubMesh\n",
      "negative particle is of type Uniform1DSubMesh\n",
      "positive particle is of type Uniform1DSubMesh\n",
      "negative particle size is of type Uniform1DSubMesh\n",
      "positive particle size is of type Uniform1DSubMesh\n",
      "current collector is of type SubMesh0D\n",
      "x_n has 20 mesh points\n",
      "x_s has 20 mesh points\n",
      "x_p has 20 mesh points\n",
      "r_n has 20 mesh points\n",
      "r_p has 20 mesh points\n",
      "y has 10 mesh points\n",
      "z has 10 mesh points\n",
      "R_n has 30 mesh points\n",
      "R_p has 30 mesh points\n"
     ]
    }
   ],
   "source": [
    "for k, t in model.default_submesh_types.items():\n",
    "    print(k, \"is of type\", t.__name__)\n",
    "for var, npts in model.default_var_pts.items():\n",
    "    print(var, \"has\", npts, \"mesh points\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solve\n",
    "Now solve the MPM with the default parameters and size distributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f79e30dbc81445a09cc42f3d8f6565fd",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, description='t', max=3581.95943869543, step=35.8195943869543), Ou…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<pybamm.plotting.quick_plot.QuickPlot at 0x7fd39283a7f0>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim = pybamm.Simulation(model)\n",
    "sim.solve(t_eval=[0, 3600])\n",
    "\n",
    "# plot some variables that depend on R\n",
    "output_variables = [\n",
    "    \"X-averaged negative particle surface concentration distribution [mol.m-3]\",\n",
    "    \"X-averaged positive particle surface concentration distribution [mol.m-3]\",\n",
    "    \"X-averaged positive electrode interfacial current density distribution [A.m-2]\",\n",
    "    \"X-averaged negative area-weighted particle-size distribution [m-1]\",\n",
    "    \"X-averaged positive area-weighted particle-size distribution [m-1]\",\n",
    "    \"Voltage [V]\",\n",
    "]\n",
    "\n",
    "sim.plot(output_variables=output_variables)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also visualise the concentration within the particles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmQAAADkCAYAAAA2PjhNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAvIUlEQVR4nO3de7xcZX3v8c+XJNzDNYDcUUQRFaggYNWKghJoeaFWj4AFpVJKT2nrafWAtrUobcXSVulLODQHkdoKVBEkehDwBqiAEjQQQsCGgBCCQrjIRQsk+3f+eJ5J1p7MWrP2zN4zs7K/b17rxcy6Pmut2b8861nPRRGBmZmZmQ3PBsNOgJmZmdl05wyZmZmZ2ZA5Q2ZmZmY2ZM6QmZmZmQ2ZM2RmZmZmQ+YMmZmZmdmQOUPWUJJ2k/SMpBlTsO9nJL1ksvc7SJLeKOmeYaejF5IukPTXNda7XtLJg0iT2VSZTrFM0mJJh1Ys/4ak9w0uRZOj7j2UdKik5YNKV9M0MkMm6XhJC/IP4OH8I37DsNPViaQ9JIWkmX3u535Jh7e+R8QDEbF5RKzuP5Xj5f0um+z9TqV8jV/a+h4R34uIlw8zTXVIer+k7xfnRcSpEXHWsNJkg+NYNr1iWUS8MiKuB5B0pqT/aFt+ZET821ASNwGDvIfTSeMyZJL+HPgM8PfADsBuwPnAMUNMVl/6DXDWTL7v05tjmTWN7+8Ui4jGTMCWwDPAuyvW2YgU5Fbk6TPARnnZocBy4C+AR4CHgZMK224C/BPwM+CXwPeBTfKyQ4CbgCeB24FDC9tdD5wF/AB4GrgOmJOXPQBETvczwOuA9+d1Pw08DvwtsCfwHeAxYCXwRWCrvI9/B8aAX+d9/G9gj7zfmXmdnYD5eX9LgT8opO9M4EvAF3L6FgMHVlzDAF6aP18MnAf8v7ztD4E9S7Zrpel9+bxXAn9ZWL4BcAZwbz7PLwHbFJafmK/9Y8BfA/cDh+dlBwE35+v/MPBZYMO87MZ83Gfz9XlP617n5WcAl7el9VzgXwq/q8/l/T6U78eMknM8E7gc+M98PX4M7FdY3jq/p4G7gHcUlrXf968A/w2szul+snDN/7aw3THAQuCpvO+5hd/dyYX1fh9YAjwBXAvsnucrH/MR0u/6DuBVw/57ns4TjmVNiWWn5Gv/MPAXNe/NHODr+fo+DnwP2CAvux84HJgLPA+8kK/D7YXrf3Le/5MU/k6B7fJ12z5//x1SXHgy3899u1yHPwWW5XtyTiFNpferkObTSXHjOeDSGvdwG+Dz+do8AXy1+Lst7HsnUhx8FLgP+NPCsoOABaS49wvgn4f9dzvlcWHYCZhQYtOPeFXrppes8wngFmD7/AO+CTir8GNYldeZBRwF/ArYOi8/L/9B7AzMAH4z/2HsnH+sR5EyFW/N37cr/BHdC7yMFAivB87Oy8b9UPO89+d0/AkwM2/z0rzfjXK6bwQ+0/ZHcXjhe/sfwA2kp+uNgf3zD/ywvOxM0j/8R+Xz+iRwS8U1bA9ij+c/jpmkP9bLSrZrpen/5nPaj/QH/Iq8/IP53uySz/NfgUvzsn1If9xvADYE/pEUrFoZsgNI/5DMzMdZAnywU5oL97qVIds93+ct8vcZpAB7SP7+1ZyWzUi/mx8Bf1hyjmfmdL2L9Bv6ECmQzMrL300KMhuQMobPAjtW3Pf3A99vO8bF5AxZvu6/JP02NiD9Fvcu/O5Ozp/fTvrH6xV5338F3JSXHQHcBmxFypy9opUmT45lOJZVxbJLSXHh1Tkdh9e4N58ELsj3ZRbwRkDt557P5T/ajns9a/+mLwL+rrDsj4Fr8ufXkDLiB+fr8L68740qrsN3SRml3YCfFo5T534tBHZlbaa+2z38f6SH1q3zNXhT4XfbissbkOLSx0gx/yWkDOMRefnNwAn58+bkeL0+T0NPwIQSC+8Fft5lnXuBowrfjwDuL/wYfs34gPII6R/6DfKy/Trs83Tg39vmXQu8L3++HvirwrL/WfjDGfdDzfPeDzzQ5TzeDvyk8L30DyD/oawGZheWfxK4OH8+E/hWYdk+wK8rjt0exC4sLDsKuLtku1aadinM+xFwbP68hBxY8/cdSZmbmfmP8tLCsk1JT5CHlxzrg8CVndJcuNfFJ7HvAyfmz28F7s2fdyBlGjcprHsc8N2S455J4R+A/Lt5GHhjyfoLgWPK7jvdM2T/Cny6ZN/XszaofgP4QFu6fkXKjL6FFIAPIT8VexruhGNZU2LZ3oV5/wB8rsa9+QRwFYV41Onc6Z4hOxxYVlj2A9bGsP9DzgAWlt9DzviUXIe5bff12xO4X79fdh4d7uGOpBK0rTvs+1DWZsgObv/tAB8BPp8/3wh8nFxCOx2mptUhewyY0+U99k6kYvqWn+V5a/YREasK339Fyn3PIT2R3dthn7sD75b0ZGsileTsWFjn5x32WeXB4hdJ20u6TNJDkp4C/iOnqY6dgMcj4unCvJ+RnobL0rfxBOoDTPTcytbfHbiycA2XkILvDvkc1lyTiPgV6X4DIOllkr4u6ef5+vw99a8PwCWkjBbA8fl7K02zgIcL6fpX0pNvmWI6x0ivjnbK6TxR0sLCvl7Vls5x972GXen8m2y3O3Bu4biPk0rDdo6I75Be8Z4H/ELSPElbTDAdNrkcyzobtVhWPLfi9a+6N+eQSquvk7RM0hk109buO8Amkg6WtDuptPDKvGx34C/a7uOujP991DqXmvdrInFrV9I9fKLLersDO7Wdw0dJ/x4AfIBUUnu3pFsl/c4E0tBITcuQ3Uwqrn57xTorSDe6Zbc8r5uVed97dlj2IOmpcqvCtFlEnF1jv1Fz/ifzvH0jYgvg90j/oHbbD6Tz20bS7MK83Uj1oUbJg8CRbddx44h4iFTKtEtrRUmbANsWtv0/wN3AXvn6fJTx16ebLwOHStoFeAdrM2QPkkrI5hTStEVEvLJiX7sW0rlBTveKHDT/L3AasG1EbAXcSfV9rLqvrfR1+k12Wu8P267tJhFxE0BE/EtEHAC8khTkPlxjnzZ1HMs6G7VYtmvhc/H6l96biHg6Iv4iIl4CHA38uaTDOuy78m8/P+x9ifQgeTzw9UJG9UHS68zifdw0Ii7t4Vy63a9Oaa1K+4Oke7hVxTqt9e5rO4fZEXEUQET8V0QcR3o4/hRwuaTNuuyz0RqVIYuIX5JebZ0n6e2SNpU0S9KRkv4hr3Yp8FeStpM0J6//H2X7LOx7jPTO/p8l7SRphqTXSdoob3+0pCPy/I1zfyq7VO8VSPUOxkjvx6vMJlfslrQz6/6D+YuyfUTEg6Q6DJ/MaduX9HTxxRrpG6QLgL/LGRfyPTomL7ucdI1/U9KGpKLqYlCYTarc+YykvYE/att36fUBiIhHSa8DPk8KAkvy/IdJFZf/SdIWkjaQtKekN1WcxwGS3pmfyj9IytDdQqprEqR7jqSTSCVkVX4B7JLPuZPPASdJOiynbed8/u0uAD4i6ZX52FtKenf+/Nr8lD2LVKet1ZDAhsSxrDGx7K/zvXklcBKpXhRU3BtJvyPppZJEilmr6fz39gtgj/xQV+YSUl3U97L2IRLSg9+p+e9akjaT9NttGdl2H5a0taRdgT8rnEu3+9VJ1T18mFSF4vx8vFmSfqvDqj8CnpJ0uqRN8u/xVZJeCyDp9yRtl3/PT+Zt1uu41agMGUBE/DPw56RKy4+SctmnkSpmQ2rls4DUImQRqRXc39bc/YfyNreSXvl8ilTn5kFSS7ePFo75YWpcv/zq7e+AH+Ri2UNKVv04qaLmL0kVIq9oW/5JUgB4UtKHOmx/HOk9/gpSsfbfRMQ3u6VvwM4ltZ66TtLTpEzMwQARsZhUMfgyUmnZ06Q6Mc/lbT9Eekp8mhSM/pPxzgT+LV+f/1Fy/EtI9TIuaZt/IqlS6V2kFkGXM/4VTrurSEHyCeAE4J0R8UJE3EVq2XYzKWC9mlTvo8p3SC3Ffi5pZfvCiPgR6R+CT5N+Gzcw/sm8td6VpN/rZfm1w53AkXnxFqRr9gRrW7H+Y5d02RRzLGtELLuB9Prx28A/RsR1eX7VvdkL+BYpk3MzcH7kvsfafDn//zFJP+508Ij4IekhaidSJqc1fwHwB6SqCE/kNL6/y7lcRapEv5B0Xz6X53e7X510u4cnkOoH302K4x/scG6rSSWI+5MaRq0ELiS1QIbU8GWxpGdI/3YcGxH/XSNtjdVq+WE2UiRtTnoq2isi7htyctaQdCapsu7vDTstZjY1JO3B2tbTq7qsPvIkBSmWLh12Wqxc40rIbP0l6ej8emAzUgnOIlJrHjMzs/WaM2Q2So5hbUeLe5GKqF2E22CS5kq6R9JSdWhtluu6fU3S7Urj/J00jHSamXVSI4ZtLelKSXdI+pGkVxWW3S9pkVLL+wVdj+V/78xsKigNNPxTUr9vy0n1mY7Lde1a63wU2DIiTpe0HakvpRdFxPPDSLOZWUvNGHYO8ExEfDw3uDovIg7Ly+4njSSxTv3gTlxCZmZT5SBgaUQsyxmsy1h3nMYAZucWaZuTKqA3vs6Oma0X6sSwfUiNPoiIu0ktZ3egB40fKHTONjNij11nDTsZZiPvtjueWxkR201kmyPevFk89nh5S/Pb7nhuMakbjZZ5ETEvf96Z8R1KLie3qi34LKnl7QpS8/v35Gbu08KG2ig2Zr3uWslsUjzNE4OOX1Avht0OvBP4vqSDSK3gdyG1tA9SrwIB/GvbvtfR+AzZHrvO4kfX7jbsZJiNvBk7/tfPuq813srHV/PDa8u7qJq1473/HREHlizu1HFvex2JI0jN8N9C6sj0m5K+FxFPTTStTbQxm3Fwxz5DzazoW3H5oOMX1IthZ5NGSVlIaoj2E9aW8r8+IlZI2p4U2+6OiBvLDtb4DJmZTZ0geCF67otxOeN7B9+FdXuaP4k0eHUASyXdB+xN6jTSzKxnfcYvqBHD8sPjSQC56sV9eSIiWqM3PCLpStIr0NIMmeuQmVmlsYr/urgV2EvSi5VGIjiW9Hqy6AGgVQF2B+DlwLJJPgUzm6b6iF9QI4ZJ2kprR1o5GbgxIp7KoyfMzutsBryN1GF3KZeQmVmp9ITZW5WuiFgl6TTgWmAGcFFELJZ0al5+AXAWcLGkRaTXA6fXbZFkZlaln/gFtWPYK4AvSFpNGu3lA3nzHYArU6EZM4FLIuKaquM5Q2ZmpQJ4od6TZOftI64Grm6bd0Hh8wrSk6OZ2aTqN35BrRh2M6nfzPbtlgH7TeRYzpCZWakAVruvQjNroKbFL2fIzKxUELywTqMiM7PR17T45QyZmZULWN2ceGZmtlbD4pczZGZWKhAvdOyKx8xstDUtfjlDZmalAhhr0BOmmVlL0+KXM2RmViqA591doZk1UNPilzNkZlYqgBeiOQHNzKylafHLGTIzKxWI1Q16wjQza2la/HKGzMxKNe0J08yspWnxyxkyM6sgVjcooJmZrdWs+OUMmZmVSkOPzBh2MszMJqxp8csZMjMrFdGsJ0wzs5amxS9nyMysVNOeMM3MWpoWv5whM7NSgXghHCbMrHmaFr+aU5ZnZkOxOlQ6mZmNsn7jl6S5ku6RtFTSGR2Wby3pSkl3SPqRpFfV3bZdc7KOZjZwTXvCNDNr6Td+SZoBnAe8FVgO3CppfkTcVVjto8DCiHiHpL3z+ofV3HYcl5CZWakAVrNB6dRNjafLD0tamKc7Ja2WtM1UnIuZTS/9xi/gIGBpRCyLiOeBy4Bj2tbZB/g2QETcDewhaYea247jDJmZlUpPmDNKpyqFJ8QjSUHrOEn7jNt/xDkRsX9E7A98BLghIh6fmrMxs+mkn/iV7Qw8WPi+PM8ruh14J4Ckg4DdgV1qbjuO30WYWakI+mk2vuYJEUBS6wmxrMj+OODSXg9mZlZUI37NkbSg8H1eRMwrfO9U0Szavp8NnCtpIbAI+Amwqua24zhDZmalWk+YFaoCWqcnxIM77UTSpsBc4LQ+kmtmtkaN+LUyIg6sWL4c2LXwfRdgxbhjRDwFnAQgScB9edq027btnCEzs1J9BrSJPCEeDfzAryvNbLLUiF/d3ArsJenFwEPAscDxxRUkbQX8KtcTOxm4MSKektR123YDrUNWo4LvlpK+Jul2SYslnTTI9JnZuvqoFNv16bLgWEb8daXjl1nz9FOpPyJWkUrtrwWWAF+KiMWSTpV0al7tFcBiSXeT6sv+WdW2VccbWAlZzSagfwzcFRFHS9oOuEfSF3PO08wGrM8nzFpPiJK2BN4E/F6vB5pqjl9mzTMJJWRExNXA1W3zLih8vhnYq+62VQb5yrJOBd8AZuf3sJsDj5Mqx5nZEAQw1mOl/ohYJan1hDgDuKj1dJmXt4LaO4DrIuLZSUjyVHH8MmuYfuLXMAwyQ1angu9ngfmk1xqzgfdExFj7jiSdApwCsNvOrgZnNlUi+nvC7PZ0mb9fDFzc80EGY0ri18ZsOiWJNbP+49egDTLrWKeC7xHAQmAnYH/gs5K2WGejiHkRcWBEHLjdts252GZNtDo2KJ2mkSmJX7PYaLLTaWYFTYpfg0xRnQq+JwFXRLKU1HR07wGlz8zaTELHiusLxy+zhmla/BpkhmxNBV9JG5Iq+M5vW+cB4DCAPPTAy4FlA0yjmRUENCqgTSHHL7OGaVr8GlgFrJoVfM8CLpa0iPSK4PSIWDmoNJrZeIEYi05v66YXxy+z5mla/BpojfgazUdXAG8bZJrMrFwEI/kkOQyOX2bN0rT45SaKZlapSU+YZmZFTYpfzpCZWanJ6FjRzGwYmha/nCEzs1KpY8XmPGGambU0LX45Q2ZmFcSqBj1hmpmt1az45QyZmZWKgBfGmhPQzMxamha/nCEzs1JNazZuZtbStPjlDJmZlQpg1QgOMWJm1k3T4pczZGZWaaxBAc3MrKhJ8as5KTWzgYsQq2KD0snMbFRNRvySNFfSPZKWSjqjw/ItJX1N0u2SFks6qbDsfkmLJC2UtKDbsVxCZmaVmlQHw8ysqJ/4JWkGcB7wVmA5cKuk+RFxV2G1PwbuioijJW0H3CPpixHxfF7+5rpDqDlDZmalAlg15pIwM2ueSYhfBwFLI2IZgKTLgGOAYoYsgNmSBGwOPA6s6uVgjrRmViror8i/W3F/XufQXKS/WNINk34SZjYt1YhfcyQtKEyntO1iZ+DBwvfleV7RZ4FXACuARcCfRcTYmiTAdZJu67DvdbiEzMzKRe9F/nWK+yVtBZwPzI2IByRt33+izcyoE79WRsSBFcs7bRxt348AFgJvAfYEvinpexHxFPD6iFiR49o3Jd0dETeWHcwlZGZWqlXkXzZ1saa4P9enaBX3Fx0PXBERDwBExCOTfQ5mNj31Gb8gPUjuWvi+C6kkrOgkUgyLiFgK3AfsDRARK/L/HwGuJMXEUs6QmVmpVseKZRPVRf51ivtfBmwt6fpcrH/i1J6RmU0XNeJXN7cCe0l6saQNgWOB+W3rPAAcBiBpB+DlwDJJm0manedvBrwNuLPqYH5laWaVVlfXFasq8q9T3D8TOIAU0DYBbpZ0S0T8dMIJNTNr0yV+VYqIVZJOA64FZgAXRcRiSafm5RcAZwEXS1pEinmnR8RKSS8Brkx1/ZkJXBIR11QdzxkyMysVfdQho15x/3JSpu5Z4FlJNwL7Ac6QmVlf+oxfeR9xNXB127wLCp9XkEq/2rdbRopltfmVpZlVEKvHNiiduqhT3H8V8EZJMyVtChwMLJn00zCzaaiv+DVwLiEzs1IBPQeuOsX9EbFE0jXAHcAYcGFEVNazMDOro5/4NQzOkJlZuUjF/j1v3qW4P38/Bzin96OYmXXQZ/waNGfIzKxU0F+lWDOzYWla/HKGzMwq1G4ebmY2YpoVv5whM7NKY2PNCWhmZkVNil/OkJlZqQiIBj1hmpm1NC1+OUNmZpVWN+gJ08ysqEnxyxkyMysViLEGNRs3M2tpWvwaaEolzZV0j6Slks4oWedQSQslLZZ0wyDTZ2brioppOnH8MmueJsWvgZWQSZoBnAe8lTRcyq2S5kfEXYV1tgLOB+ZGxAOSth9U+sysg4BoUJH/VHH8MmughsWvQZaQHQQsjYhlEfE8cBlwTNs6xwNXRMQDABHxyADTZ2YdRKh0mkYcv8waqEnxa5AZsp2BBwvfl+d5RS8DtpZ0vaTbJJ3YaUeSTpG0QNKCRx9bPUXJNbMgNRsvm6aRKYlfL/DcFCXXzJoWvwZZqb/T2be/xp0JHAAcBmwC3Czploj46biNIuYB8wAO3G/jUXwVbLZ+CGAEnySHYEri1xbaxvHLbKo0LH4NsoRsObBr4fsuwIoO61wTEc9GxErgRmC/AaXPzDqIsfJpGnH8MmugfuNXt8Y8kraU9DVJt+fGPCfV3bbdIDNktwJ7SXqxpA2BY4H5betcBbxR0kxJmwIHA0sGmEYzG0fEWPk0jTh+mTVOf/Gr0JjnSGAf4DhJ+7St9sfAXRGxH3Ao8E+SNqy57TgDe2UZEasknQZcC8wALoqIxZJOzcsviIglkq4B7gDGgAsj4s5BpdHM2jSsp+up4vhl1kD9x681jXkAJLUa89xVWCeA2ZIEbA48DqwiPZB123acgXYMGxFXA1e3zbug7fs5wDmDTJeZVXCGDHD8Mmuk6vg1R9KCwvd5uY5nS6fGPAe37eOzpNLyFcBs4D0RMSapzrbjuKd+M6u2nlU7l7RtRDw27HSY2QBUx6+VEXFgxfI6jXmOABYCbwH2BL4p6Xs1tx2npwxZLnp7IX99OCL+dy/7MbMRF0AfdcUkzQXOJb3muzAizm5bfiip7tV9edYVEfGJng/YPT1fAR6RtEVOz3en6lhmNmR9xi/qNeY5CTg7IgJYKuk+YO+a247TawnZzRFxLqSnzR73YWYNED2WkNXp3T77XkT8Tl+JrO/uiPjLnL7zAGfIzNZjvcavbE1jHuAhUmOe49vWeYDU1c33JO0AvBxYBjxZY9txes2QHSNpDLi2vY8dM1vP9P6EWadC7KDNlfQ4cDup4q2Zrc/6KCGr05gHOAu4WNIi0mvK03O3N3Tatup4vWbITiD1r/O7kvaMiJN73I+ZjbIA9d7fWN1Kra+TdDupOP9D3YJWn+YCrwVeD2wn6d8i4n1TeDwzG5b+4lfaRZfGPBGxAnhb3W2r9JQhi4iHSEVwtQ9kZk2kflop1anU+mNg94h4RtJRwFeBvXpNbTcR8Sgpbjl2ma33usavkdK1Y1hJJ0h6VNJySe/L8w6R9LeSbpv6JJrZUI1VTLmVUmEqNhnvWqk1Ip6KiGfy56uBWZLmTEayHbvMrEv8Gil1eur/GHAUsD/wYknfBL4MbAh8cMpSZmajISqmal17t5f0otyhIpIOIsWkyeqSwrHLbLrrPX4NXJ1Xls9ExK0Akj4O/AJ4WUQ8OZUJM7MREKAeK8XWrBD7LuCPJK0Cfg0cm5uPTwbHLrPprI/4NQx1MmQvknQKcE+eljugmU0jfWSPalSI/Sypp+up4NhlNt2NYElYmToZsr8B9gXeC7yaNGbTt4CfAD+JiEumMH1mNmRqUEBr49hlNs01KX51zZC1VdJF0i6kIPdq0ijmDmpm66v+e7oeGscus2muYfGra4ZM0uuAW1r1OiJiOan1lJuNm00HDXrCLHLsMrMmxa86rSzfB/xY0mWS3i/pRVOdKDMbHRorn0acY5fZNNek+FXnleWpAJL2JhXzXyxpS9IYcNcAP4iI1VOaSjMbngY9YRY5dplZk+JXnRIyACLi7oj4dETMBd4CfB94N/DDqUqcmQ2XcrPxsqkJOsSuu3HsMlvvNS1+9TqW5edJA/MG6WnTzNZXDXrCrOHfgUeBLYEPDzktZjbVGhS/es2Q3RwR5wJI2nYS02NmI2YU61r04Z6I+EsASefhB0qz9VqT4levGbJjJI0B10bETyczQWY2QqJZ/fjUMFfS48DtpFJ+M1tfNSx+1apDJmkDSR8tzDoBuBf4XUkXTknKzGw0NGhw3nYdYtdcYAnwemA7Sf82nJSZ2UD0Gb8kzZV0j6Slks7osPzDkhbm6U5JqyVtk5fdL2lRXrag27FqlZBFxJikw4G/z98fAh7C/fmYrfea9ITZrkPsepQUtxy7zKaBfuKXpBnAecBbSX0Y3ippfkTc1VonIs4BzsnrHw38r4h4vLCbN0fEyjrHq93KEviJpL+RNJFtzKzpomJqBscus+mqv/h1ELA0IpZFxPPAZcAxFesfB1zaa1InEqB2BY4FVki6StJZkt7d64HNrAGiWR0rlnDsMpuO+o9fOwMPFr4vz/PWIWlTUpWIr4xPAddJuk3SKd0OVrtSf0T8j3zQjYBXksaDOwj4ct19mFkDNackrCPHLrNprDp+zWmr2zWvbQzcTp2Vle3xaFJn08XXla+PiBWStge+KenuiLixLDETbmUZEc8BP86Tma3HRKNKwio5dplNLzXi18qIOLBi+XJSCXvLLsCKknWPpe11ZUSsyP9/RNKVpAfB0gzZQOtUdGutUFjvtbmlwrsGmT4za5ObjZdN04njl1nD9B+/bgX2kvRiSRuSMl3z21fKQ7K9CbiqMG8zSbNbn4G3AXdWHazXfsgmrE5rhcJ6nwKuHVTazKzCelJC1g/HL7OG6iN+RcQqSaeR/p5nABdFxGJJp+blF+RV3wFcFxHPFjbfAbhSEqS81iURcU3V8QaWIaPQWgFAUqu1wl1t6/0JqVLcaweYNjMr0c8rS0lzgXNJwezCiDi7ZL3XArcA74mIy3s/4pRx/DJroH6rXETEOt3kFDJire8XAxe3zVsG7DeRYw3ylWXX1gqSdiblNMedbDtJp0haIGnBo4+tnvSEmllW1WS8S5F/oVTpSGAf4DhJ+5SsN+qlSlMSv17guUlPqJllfcSvYRhkhqxOa4XPAKdHRGUuKyLmRcSBEXHgdtvOmKz0mVkHfTQbr9uHT6tU6ZFJTfjkmpL4NYuNJit9ZtZBk7rtGeQryzqtFQ4ELsvvXOcAR0laFRFfHUgKzWwdXSq/VjUb71SqdPC4fa8tVXoLo/2az/HLrIGa1PhokBmyNa0VSMMuHQscX1whIl7c+izpYuDrDmZmQxR0qxRb1Wx8QqVKOSMzqhy/zJqme/waKQPLkE2gtYKZjQjROVdV03pTquT4ZdY8fcavgRtkCVmt1gqF+e8fRJrMrFofdS3Wq1Ilxy+z5hnFumJlBpohM7MG6jGguVTJzIbOGTIzWy/02SO/S5XMbGgaNqKIM2RmVqlJRf5mZkVNil/OkJlZtQY9YZqZjdOg+OUMmZmVi2Y9YZqZrdGw+OUMmZlVa9ATppnZOA2KX86QmVkp0awnTDOzlqbFL2fIzKxcgMYa9IhpZtbSsPjlDJmZVWpSs3Ezs6Imxa8Nhp0AMxttGiufzMxGWb/xS9JcSfdIWirpjA7LPyxpYZ7ulLRa0jZ1tm3nDJmZVYuKycxslPURvyTNAM4DjgT2AY6TtM+43UecExH7R8T+wEeAGyLi8TrbtnOGzMzKhUvIzKyh+o9fBwFLI2JZRDwPXAYcU7H+ccClPW7rDJmZlROpDkbZZGY2qmrErzmSFhSmU9p2sTPwYOH78jxv3WNJmwJzga9MdNsWV+o3s0pNaqVkZlbUJX6tjIgDqzbvMK9sh0cDP4iIx3vYFnCGzMyqBGj1sBNhZtaD/uPXcmDXwvddgBUl6x7L2teVE90W8CtLM+vGlfrNrKn6i1+3AntJerGkDUmZrvntK0naEngTcNVEty1yCZmZlWtYx4pmZmv0Gb8iYpWk04BrgRnARRGxWNKpefkFedV3ANdFxLPdtq06njNkZlbJlffNrKn6jV8RcTVwddu8C9q+XwxcXGfbKn5laWalWmPB9dpsvEanisdIuiN3qrhA0hum4DTMbBrqN34NmkvIzKxcRJp6UOgY8a2kCq63SpofEXcVVvs2MD8iQtK+wJeAvftMtZlZX/FrGJwhM7NKfTxJrukYEUBSq2PENRmyiHimsP5muKmAmU2iUSwJK+MMmZmVC9DqyjzSHEkLCt/nRcS8/LlTx4gHt+9A0juATwLbA7/dX4LNzLLu8WukOENmZtWq41lVx4q1OkaMiCuBKyX9FnAWcPhEk2hm1lFz8mPOkJlZtT6ajU+oY8SIuFHSnpLmRMTKXg9qZtbSpG573MrSzCr1MZZl144RJb1UkvLn1wAbAo9N/lmY2XTUpLF4B5ohq9EE/r25Cfwdkm6StN8g02dm4yl3rFg2VYmIVUCrY8QlwJdanSq2OlYEfhe4U9JCUovM90SMZrMoxy+zZuknfg3DwF5Z1mwCfx/wpoh4QtKRwDw6VAI2swHqo5VSt04VI+JTwKd6P8JgOH6ZNZRbWXZUpwn8TYX1byHVOTGzYfHQSS2OX2ZN07D4NchXlp2awO9csf4HgG90WiDplNyr94JHH+tvKHczq1Je3N+kQDcJpiR+vcBzk5hEMxuvWfFrkCVktZrAA0h6MymgdRxGJfdzNA/gwP02Hr2rarY+Gc0qXYM2JfFrC23ji2s2lRoUvwaZIavVBD4Pn3IhcGREuLWV2TA1rGPFKeT4ZdY0DYtfg3xlWacJ/G7AFcAJEfHTAabNzMpExTR9OH6ZNVGD4tfASsgiYpWkVhP4GcBFrSbwefkFwMeAbYHzc9dEqyp6ATezAdBYg5opTRHHL7Nm6jd+SZoLnEv6u78wIs7usM6hwGeAWaTRS96U598PPA2spkY8GGhP/TWawJ8MnDzINJlZhaBRzcankuOXWcP0Gb/qdHcjaSvgfGBuRDwgafu23by57sgjHjrJzEqJcAmZmTXSJMSvrt3dAMcDV0TEAwAR8UivB/PQSWZWLoDVUT6ZmY2q7vFrTqsLmjyd0raHOt3dvAzYWtL1km6TdGJbCq7L89v3vQ6XkJlZJTWo2biZWVGX+LWyS72uOt3dzAQOAA4DNgFulnRLbtjz+ohYkV9jflPS3RFxY9nBnCEzswoBfmVpZo3Ud/yq093NclLG7lngWUk3AvsBP42IFZBeY0q6kvQKtDRD5leWZlYuSB0rlk1mZqOq//jVtbsb4CrgjZJmStqUNH7tEkmbSZoNIGkz4G3AnVUHcwmZmVVqUseKZmZF/cSvOt3dRMQSSdcAd5DadF4YEXdKeglwZe4CZyZwSURcU3U8Z8jMrJpLwsysqfqMX926u8nfzwHOaZu3jPTqsjZnyMysXASsdh0yM2ughsUvZ8jMrJor9ZtZUzUofrlSv5mVC2AsyqcuJM2VdI+kpZLO6LD8vZLuyNNNkiZUxG9mVqrP+DVoLiEzswoBY6t72rLOsCPAfcCbIuIJSUcC80itlMzM+tR7/BoGZ8jMrFzrCbM3XYcdiYibCuvfQurnx8ysf/3Fr4FzhszMqlXXwZgjaUHh+7yImJc/dxp2pKr06wPAN3pKo5lZJw2qQ+YMmZlV6NqBYtXQI3WGHUkrSm8mZcjeMLH0mZmVaVYH1s6QmVm5AFb3XAejzrAjSNoXuBA4MiIe6/VgZmbj9Be/Bs4ZMjOr0Fc/PmuGHQEeIg07cnxxBUm7AVcAJ+TBeM3MJon7ITOz9UVARG8Brc6wI8DHgG2B8/MQI6sqXoGamdXXR/waBmfIzKxaH0+Y3YYdiYiTgZN7PoCZWRWXkJnZeiGiUa2UzMzWaFj8cobMzCpFgyrFmpkVNSl+OUNmZhWa1WzczGytZsUvj2VpZuVazcbLJjOzUTUJ8avbeLx5nUMlLZS0WNINE9m2yCVkZlYqIhpV5G9m1tJv/KozHq+krYDzgbkR8YCk7etu284lZGZWKcaidDIzG2V9xq814/FGxPNAazzeouOBKyLiAYCIeGQC247T+BKy2+54buWMHf/rZ8NORwdzgJXDTkQNTufkGuV07j7RDZ7miWu/NfalORWrjOq5NsLTPLHyW3H5KMYvGO3fcpHTOblGNZ1TEb82rhiLF+qNx/syYJak64HZwLkR8YWa247T+AxZRGw37DR0ImlBEzq4dDonV1PSWVdEzB12GtZnoxq/oDm/ZadzcjUlnXVMQvyqMx7vTOAA4DBgE+BmSbfU3HadHZmZmZnZeHXG410OrIyIZ4FnJd0I7Fdz23Fch8zMzMxsXWvG45W0IWk83vlt61wFvFHSTEmbkl5LLqm57TguIZs687qvMhKczsnVlHSaddOU37LTObmaks4pV2c83ohYIuka4A5gDLgwIu4E6LRt1fEUDeo0zczMzGx95FeWZmZmZkPmDJmZmZnZkDlD1qduQyPkIRV+mYdVWCjpY0NK50WSHpF0Z8lySfqXfB53SHrNoNOY09EtnUO/npJ2lfRdSUvyUBl/1mGdkbieZlUcvyZXE+JXTodj2CiKCE89TqSKevcCLwE2BG4H9mlb51Dg6yOQ1t8CXgPcWbL8KOAbpL5TDgF+OKLpHPr1BHYEXpM/zwZ+2uG+j8T19OSpbHL8Gko6R+V6OoaN4OQSsv5MeGiEYYmIG4HHK1Y5BvhCJLcAW0nacTCpW6tGOocuIh6OiB/nz0+Tmjjv3LbaSFxPswqOX5OsCfELHMNGlTNk/ek0NEL7jxrgdZJul/QNSa8cTNImrO65jIKRuZ6S9gB+A/hh26ImXU+bnhy/hmOkrqdj2OhwP2T9qTM0wo+B3SPiGUlHAV8F9prqhPVgwsM8DMnIXE9JmwNfAT4YEU+1L+6wySheT5u+HL8Gb6Sup2PYaHEJWX+6Do0QEU9FxDP589WkQUirBjsdlgkP8zAMo3I9Jc0iBbIvRsQVHVZpxPW0ac3xa8BG6Xo6ho0eZ8j603VoBEkvkqT8+SDSNX9s4Cntbj5wYm5Zcwjwy4h4eNiJajcK1zMf/3PAkoj455LVGnE9bVpz/BqwUbmejmGjya8s+xA1hlUA3gX8kaRVwK+BYyNi4MW+ki4ltfCZI2k58DfArEI6rya1qlkK/Ao4adBprJnOUbierwdOABZJWpjnfRTYrZDOkbieZmUcv4aSzpG4njiGjSQPnWRmZmY2ZH5laWZmZjZkzpCZmZmZDZkzZGZmZmZD5gyZmZmZ2ZA5Q2Y2ZN0GJO5hf6sLgxfP776FmVlvHL8mj1tZmg2ZpN8CniGNG/eqSdjfMxGxef8pMzOr5vg1eVxCZjZknQYklrSnpGsk3Sbpe5L2HlLyzMxKOX5NHmfIbA1Jfyjp53ng23slnVhjmz0k/brQuWC/adgkF1U/P6JDtAzKPOBPIuIA4EPA+RPYdmNJCyTdIuntU5I6sxHkGDYyHL964J76rWhf4MyIuCAP63E18IUa290bEftPRgIi4tfA/pLun4z9NZHSgL+/CXw5j7ICsFFe9k7gEx02eygijsifd4uIFZJeAnxH0qKIuHeq0202AhzDhszxq3fOkFnRq4HL8+f7gOd72Ymk64E/jIh7JG0L3BARr5K0B3AN8H3gEOB24PPAx4HtgfdGxI/6OoP1wwbAk53+gciDAHcaCLi4zor8/2X5XvwGMC0Cmk17jmHD5/jVI7+ytKJXA/fkgWdPA/6yx/28FPiv/HlfYFHbsnPz/L2B44E3kIq1P9rj8dYrEfEUcJ+kd0MaCFjSfnW2lbS1pNbT6BzSmHV3TVlizUaLY9iQOX71zhkyA0DSrsBsUhH/I6QAc3EP+9mdVPw8lmftC9xRWOW+iFiUly8Gvp0H110E7NHzCTRYHpD4ZuDlkpZL+gDwXuADkm4nXadjau7uFcCCvN13gbMjYtoENJu+HMOGw/Fr8viVpbXsC9wYEW+RtDVwJ/A6SQ8A/w7MBw6JiPd02c/+jA9eBwD/Wfj+XOHzWOH7GNP09xgRx5UsmtvDvm4ilRKYTTdlMWwp6TXjtcDLgXcVMlud7I9jWG2OX5PHJWTW8mrgJwAR8QRwCfDbwH7AVyPi08CqGvvZD9gYQNJepCejRZVbmJn1ryyGvRa4NCI+Qio527bLfhzDbCicIbOWNcEs+xpwFCk4XZvn1elFeH9gg1zk/DFgCfC+yUummVlHZTHstaTK9wBbRsSjXfazP45hNgTuqd8qSboIOBnYBjgjIj7UtnwP4OutHprz64HfiIin+zzu/cCBEbGyn/2Y2fSW6zitJL1OvDwivt22fA8cw2wETLv33TYxEfH7+eNKUiuidquBLXOnim8ExvoJZJI2IVUQnUWqk2Fm1o9VEfEnFcsdw2wkuITMzMzMbMhch8zMzMxsyJwhMzMzMxsyZ8jMzMzMhswZMjMzM7Mhc4bMzMzMbMicITMzMzMbMmfIzMzMzIbMGTIzMzOzIfv/jaUsw561AO0AAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 720x216 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Concentrations as a function of t, r and R\n",
    "c_s_n = sim.solution[\n",
    "    \"X-averaged negative particle concentration distribution [mol.m-3]\"\n",
    "]\n",
    "c_s_p = sim.solution[\n",
    "    \"X-averaged positive particle concentration distribution [mol.m-3]\"\n",
    "]\n",
    "\n",
    "# r_n, r_p\n",
    "r_n = sim.solution[\"r_n [m]\"].entries[:, 0, 0]\n",
    "r_p = sim.solution[\"r_p [m]\"].entries[:, 0, 0]\n",
    "# dimensional R_n, R_p\n",
    "R_n = sim.solution[\"Negative particle sizes [m]\"].entries[:, 0]\n",
    "R_p = sim.solution[\"Positive particle sizes [m]\"].entries[:, 0]\n",
    "t = sim.solution[\"Time [s]\"].entries\n",
    "\n",
    "\n",
    "def plot_concentrations(t):\n",
    "    f, axs = plt.subplots(1, 2, figsize=(10, 3))\n",
    "    plot_c_n = axs[0].pcolormesh(\n",
    "        R_n, r_n, c_s_n(r=r_n, R=R_n, t=t), vmin=0.15, vmax=0.8\n",
    "    )\n",
    "    plot_c_p = axs[1].pcolormesh(\n",
    "        R_p, r_p, c_s_p(r=r_p, R=R_p, t=t), vmin=0.6, vmax=0.95\n",
    "    )\n",
    "    axs[0].set_xlabel(r\"$R_n$ [$\\mu$m]\")\n",
    "    axs[1].set_xlabel(r\"$R_p$ [$\\mu$m]\")\n",
    "    axs[0].set_ylabel(r\"$r_n / R_n$\")\n",
    "    axs[1].set_ylabel(r\"$r_p / R_p$\")\n",
    "    axs[0].set_title(\"Concentration in negative particles [mol.m-3]\")\n",
    "    axs[1].set_title(\"Concentration in positive particles [mol.m-3]\")\n",
    "    plt.colorbar(plot_c_n, ax=axs[0])\n",
    "    plt.colorbar(plot_c_p, ax=axs[1])\n",
    "\n",
    "    plt.show()\n",
    "\n",
    "\n",
    "# initial time\n",
    "plot_concentrations(t[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmQAAADkCAYAAAA2PjhNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA2tElEQVR4nO3de7gkVXnv8e9vbtxvMkCUqyKKqEAUAaNGFJSBhAc1egS8EjmEnJATT6JHNInBmEQNiUoe4ZA5iMREIEZB0IOAUREvoAw6XIYBMwwIw6AwXIThNrP3fs8fVT1Tu+laXbv6WrN/H55+6O5Vq2p1de93Vq2q9ZYiAjMzMzMbnTmjboCZmZnZbOcOmZmZmdmIuUNmZmZmNmLukJmZmZmNmDtkZmZmZiPmDpmZmZnZiLlD1lCS9pC0VtLcAax7raTn9Xu9wyTpNZJuH3U76pB0jqS/rLDc1ZJOGkabzAZlNsUyScskHZYo/6ak9wyvRf1R9TuUdJikVcNqV9M0skMm6QRJS/IfwH35j/jVo25XJ5L2khSS5vW4nrskHdF6HRF3R8TWETHZeyuny9e7st/rHaR8Hz+/9Toivh8RLxxlm6qQ9F5JPyi+FxGnRMTHR9UmGx7HstkVyyLixRFxNYCk0yX9W1v5URHxLyNp3AwM8zucTRrXIZP0p8Bngb8DdgH2AM4Gjh1hs3rSa4CzZvL3Prs5llnT+PsdsIhozAPYDlgLvC2xzGZkQW51/vgssFledhiwCvgz4H7gPuDEQt0tgH8EfgH8GvgBsEVedijwI+AR4EbgsEK9q4GPAz8EHgOuAhbmZXcDkbd7LfBK4L35sp8BHgL+Btgb+A7wILAG+BKwfb6OfwWmgCfzdfxvYK98vfPyZZ4DXJavbwXw3wvtOx34MvDFvH3LgIMS+zCA5+fPzwfOAv5fXvfHwN4l9Vptek/+udcAf14onwOcBtyRf84vA88qlL873/cPAn8J3AUckZcdDFyb7//7gM8BC/Kya/LtPp7vn7e3vuu8/DTgK21tPRP4p8Lv6vP5eu/Nv4+5JZ/xdOArwL/n++OnwAGF8tbnewy4FXhzoaz9e/8q8BQwmbf7kcI+/5tCvWOBpcCj+boXFX53JxWW+31gOfAwcCWwZ/6+8m3eT/a7vgl4yaj/nmfzA8eypsSyk/N9fx/wZxW/m4XAN/L9+xDwfWBOXnYXcASwCFgHrM/3w42F/X9Svv5HKPydAjvl+23n/PXvksWFR/Lvc/8u++F/Aivz7+SMQptKv69Cmz9EFjeeBi6s8B0+C/hCvm8eBr5W/N0W1v0csjj4AHAn8D8LZQcDS8ji3q+AT4/673bgcWHUDZhRY7Mf8UTrSy9Z5q+B64Cd8x/wj4CPF34ME/ky84GjgSeAHfLys/I/iF2BucBv5X8Yu+Y/1qPJOhVvyF/vVPgjugN4AVkgvBr4ZF427Yeav/fevB1/DMzL6zw/X+9mebuvAT7b9kdxROF1+x/A98iOrjcHDsx/4IfnZaeT/cN/dP65PgFcl9iH7UHsofyPYx7ZH+tFJfVabfq/+Wc6gOwP+EV5+fvz72a3/HP+M3BhXrYf2R/3q4EFwD+QBatWh+zlZP+QzMu3sxx4f6c2F77rVodsz/x73jZ/PZcswB6av/5a3patyH43PwH+oOQznp63661kv6EPkAWS+Xn528iCzByyjuHjwLMT3/t7gR+0beN88g5Zvt9/TfbbmEP2W9y38Ls7KX/+JrJ/vF6Ur/svgB/lZUcCNwDbk3XOXtRqkx+OZTiWpWLZhWRx4aV5O46o8N18Ajgn/17mA68B1P7Z88/yb23bvZqNf9PnAX9bKPsj4Ir8+cvIOuKH5PvhPfm6N0vsh++SdZT2AH5e2E6V72spsDsbO/XdvsP/R3bQukO+D15b+N224vIcsrj0UbKY/zyyDuORefm1wLvy51uTx+tN+THyBsyosfAO4JddlrkDOLrw+kjgrsKP4UmmB5T7yf6hn5OXHdBhnR8C/rXtvSuB9+TPrwb+olD2Pwp/ONN+qPl77wXu7vI53gT8rPC69A8g/0OZBLYplH8COD9/fjrwn4Wy/YAnE9tuD2LnFsqOBm4rqddq026F934CHJc/X04eWPPXzybr3MzL/ygvLJRtSXYEeUTJtt4PXNKpzYXvungk9gPg3fnzNwB35M93Ies0blFY9njguyXbPZ3CPwD57+Y+4DUlyy8Fji373uneIftn4DMl676ajUH1m8D72tr1BFln9PVkAfhQ8qNiP0b7wLGsKbFs38J7fw98vsJ389fApRTiUafPTvcO2RHAykLZD9kYw/4PeQewUH47ecenZD8savtevz2D7+v3yz5Hh+/w2WQjaDt0WPdhbOyQHdL+2wE+DHwhf34N8DHyEdrZ8GjaNWQPAgu7nMd+Dtkwfcsv8vc2rCMiJgqvnyDrfS8kOyK7o8M69wTeJumR1oNsJOfZhWV+2WGdKfcUX0jaWdJFku6V9Cjwb3mbqngO8FBEPFZ47xdkR8Nl7dt8BtcDzPSzlS2/J3BJYR8uJwu+u+SfYcM+iYgnyL5vACS9QNI3JP0y3z9/R/X9A3ABWUcL4IT8datN84H7Cu36Z7Ij3zLFdk6RnTp6Tt7Od0taWljXS9raOe17r2B3Ov8m2+0JnFnY7kNko2G7RsR3yE7xngX8StJiSdvOsB3WX45lnY1bLCt+tuL+T303Z5CNVl8laaWk0yq2rd13gC0kHSJpT7LRwkvysj2BP2v7Hndn+u+j0mep+H3NJG7tTvYdPtxluT2B57R9ho+Q/XsA8D6ykdrbJF0v6Xdn0IZGalqH7Fqy4eo3JZZZTfZFt+yRv9fNmnzde3cou4fsqHL7wmOriPhkhfVGxfc/kb+3f0RsC7yT7B/UbuuB7PM9S9I2hff2ILseapzcAxzVth83j4h7yUaZdmstKGkLYMdC3f8D3Absk++fjzB9/3TzH8BhknYD3szGDtk9ZCNkCwtt2jYiXpxY1+6Fds7J2706D5r/FzgV2DEitgduIf09pr7XVvs6/SY7LfcHbft2i4j4EUBE/FNEvBx4MVmQ+2CFddrgOJZ1Nm6xbPfC8+L+L/1uIuKxiPiziHgecAzwp5IO77Du5N9+frD3ZbIDyROAbxQ6qveQnc4sfo9bRsSFNT5Lt++rU1tTbb+H7DvcPrFMa7k72z7DNhFxNEBE/FdEHE92cPwp4CuStuqyzkZrVIcsIn5NdmrrLElvkrSlpPmSjpL09/liFwJ/IWknSQvz5f+tbJ2FdU+RnbP/tKTnSJor6ZWSNsvrHyPpyPz9zfN8Krul1wpk1x1MkZ0fT9mG/MJuSbvyzH8wf1W2joi4h+wahk/kbduf7OjiSxXaN0znAH+bd1zIv6Nj87KvkO3j35K0gGyouhgUtiG7uHOtpH2BP2xbd+n+AYiIB8hOB3yBLAgsz9+/j+zC5X+UtK2kOZL2lvTaxOd4uaS35Efl7yfr0F1Hdq1JkH3nSDqRbIQs5VfAbvln7uTzwImSDs/btmv++dudA3xY0ovzbW8n6W3581fkR9nzya5pa00ksBFxLGtMLPvL/Lt5MXAi2XVRkPhuJP2upOdLElnMmqTz39uvgL3yg7oyF5Bdi/oONh5EQnbgd0r+dy1JW0n6nbaObLsPStpB0u7AnxQ+S7fvq5PUd3gf2SUUZ+fbmy/ptzss+hPgUUkfkrRF/nt8iaRXAEh6p6Sd8t/zI3mdTTpuNapDBhARnwb+lOyi5QfIetmnkl2YDdksnyVkM0JuJpsF9zcVV/+BvM71ZKd8PkV2zc09ZDPdPlLY5gepsP/yU29/C/wwH5Y9tGTRj5FdqPlrsgsiL24r/wRZAHhE0gc61D+e7Dz+arJh7b+KiG91a9+QnUk2e+oqSY+RdWIOAYiIZWQXBl9ENlr2GNk1MU/ndT9AdpT4GFkw+nemOx34l3z//LeS7V9Adl3GBW3vv5vsotJbyWYEfYXpp3DaXUoWJB8G3gW8JSLWR8StZDPbriULWC8lu+4j5TtkM8V+KWlNe2FE/ITsH4LPkP02vsf0I/PWcpeQ/V4vyk873AIclRdvS7bPHmbjLNZ/6NIuGzDHskbEsu+RnX78NvAPEXFV/n7qu9kH+E+yTs61wNmR5x5r8x/5/x+U9NNOG4+IH5MdRD2HrJPTen8J8N/JLkV4OG/je7t8lkvJLqJfSva9fD5/v9v31Um37/BdZNcH30YWx9/f4bNNko0gHkg2MWoNcC7ZDGTIJr4sk7SW7N+O4yLiqQpta6zWzA+zsSJpa7Kjon0i4s4RN2cDSaeTXaz7zlG3xcwGQ9JebJw9PdFl8bEnKchi6YpRt8XKNW6EzDZdko7JTw9sRTaCczPZbB4zM7NNmjtkNk6OZWOixX3Ihqg9hNtgkhZJul3SCnWYbZZf6/Z1STcqu8/fiaNop5lZJxVi2A6SLpF0k6SfSHpJoewuSTcrm3m/pOu2/O+dmQ2CshsN/5ws79sqsuuZjs+vtWst8xFgu4j4kKSdyHIp/UZErBtFm83MWirGsDOAtRHxsXzC1VkRcXhedhfZnSSecX1wJx4hM7NBORhYEREr8w7WRTzzPo0BbJPPSNua7AL0xl+zY2abhCoxbD+ySR9ExG1kM2d3oYbG3yh0gTaLzdUtt9/MzCS51eA3OOTWqOb2BvEZ6q6z7i6r/dkH8B0N4Ht49In71kTETjNZ3ZGv2yoefKh8pvkNNz29jCyNRsviiFicP9+V6QklV5HPqi34HNnM29Vk0+/fnk9znxUWPmtu7LX7/FE3w8ZIJFJ81T2fVXedtbeXOPMWiSCV2t6ym9cPO35BtRh2I/AW4AeSDiabBb8b2Uz7IMsqEMA/t637GRrfIdtcW3PovCP7uk7NqfuPYc0Bx9T2Ev8wa07N7aX+sU+tc255WTKVTurzpbZXtyzRzuRnnzu379uL5GcoL0ptL7nOueWf76obPvaL0sISax6a5MdXlqeomv/sO56KiINKijs1pj3mHkk2Df/1ZIlMvyXp+xHx6Ezb2kR77T6fn1y5x6ibYWNkfWJS51SiyzJF+XHM+sQxzmRinal6qaOmdYkO2fpEh2wq0SPbd4/7hh2/oFoM+yTZXVKWkk1E+xkbR/lfFRGrJe1MFttui4hryjbW+A6ZmQ1OEKyP2rkYVzE9O/huPDPT/IlkN68OYIWkO4F9yZJGmpnV1mP8ggoxLD94PBEgv/TizvxBRLTu3nC/pEvIToGWdsh8DZmZJU0l/uviemAfSc9VdieC48hOTxbdDbQugN0FeCGwss8fwcxmqR7iF1SIYZK218Y7rZwEXBMRj+Z3T9gmX2Yr4I1kCbtLeYTMzEplR5j1LumKiAlJpwJXAnOB8yJimaRT8vJzgI8D50u6mez0wIeqzkgyM0vpJX5B5Rj2IuCLkibJ7vbyvrz6LsAl2aAZ84ALIuKK1PbcITOzUgGsr3Yk2bl+xOXA5W3vnVN4vprsyNHMrK96jV9QKYZdS5Y3s73eSuCAmWzLHTIzKxXApHMVmlkDNS1+uUNmZqWCYH3tye9mZqPTtPjlDpmZlQuYbE48MzPbqGHxyx2ymRqnXGPDzieWyu81gBxezEmsc16iLJUXbF693GapepH47OntJfLxzEskT0x+vv4mqQ2UzBtkNpulc3+Vp1tI1XsqUS91gXpqJGh9olPyVJTHk/VRHtueivLuw7pkvfJEyOtJxHXuS5R11rT45Q6ZmZUK0skazczGVdPilztkZlYqgHVOV2hmDdS0+OUOmZmVCmB94pSGmdm4alr8cofMzEoFYrJBR5hmZi1Ni1/ukJlZqaYdYZqZtTQtfrlDZmYJYrJBAc3MbKNmxS93yMysVHbrkdRUdDOz8dS0+DVrO2RK5QWrqym5xmrnDKtZb17iZ5bKJ5Yoq50XbH7Ndc5PbS+RTyxRL5lrLFGWqjeVaEsdEc06wjTrt7q5xp6OifJ6yTxkqbLyv++nErm/npgqz/2VyguWyjX2+NRmiXoLyssSbXliqrweLEuUdda0+DVrO2Rm1l3TjjDNzFqaFr/cITOzUoFYnzhKNjMbV02LX80ZyzOzkZgMlT7MzMZZr/FL0iJJt0taIem0DuU7SLpE0k2SfiLpJVXrtmtO19HMhq5pR5hmZi29xi9Jc4GzgDcAq4DrJV0WEbcWFvsIsDQi3ixp33z5wyvWncYjZGZWKoBJ5pQ+uqlwdPlBSUvzxy2SJiU9axCfxcxml17jF3AwsCIiVkbEOuAi4Ni2ZfYDvg0QEbcBe0napWLdadwhM7NS2RHm3NJHSuEI8SiyoHW8pP2mrT/ijIg4MCIOBD4MfC8iHhrMpzGz2aSX+JXbFbin8HpV/l7RjcBbACQdDOwJ7Fax7jQ+F2FmpSLoZdr4hiNEAEmtI8SyIfvjgQvrbszMrKhC/FooaUnh9eKIWFx43elCs2h7/UngTElLgZuBnwETFetO4w6ZmZVqHWEmpAJapyPEQzqtRNKWwCLg1B6aa2a2QYX4tSYiDkqUrwJ2L7zeDVg9bRsRjwInAkgScGf+2LJb3XbukI2DQSR/VaJe3eSvySSuiZ9SKlFrap2pJK6JsqlkvUSi1gWpJK6JsvmpxLA1E7wm65UWMZmoV0ePAW0mR4jHAD/06UobhUEkf30qUe+pKB8oeSIxovN4IolrKqnq44lErakEr08kytL1yrf3xGR5vbWJsjoqxK9urgf2kfRc4F7gOOCE4gKStgeeyK8TOwm4JiIeldS1bruhXkNW4QLf7SR9XdKNkpZJOnGY7TOzZ+rhotiuR5cFxzHmpysdv8yap5eL+iNigmzU/kpgOfDliFgm6RRJp+SLvQhYJuk2sutl/yRVN7W9oY2QVZwC+kfArRFxjKSdgNslfSnveZrZkPV4hFnpCFHSdsBrgXfW3dCgOX6ZNU8fRsiIiMuBy9veO6fw/Fpgn6p1U4Z5yrLKBb4BbJOfh90aeIjs4jgzG4EApmpe1B8RE5JaR4hzgfNaR5d5eSuovRm4KiIe70OTB8Xxy6xheolfozDMDlmVC3w/B1xGdlpjG+DtEc+806qkk4GTATZny4E01syym/P2coTZ7egyf30+cH7tjQzHQOLXHrv6Ml6zQek1fg3bMLuOVS7wPRJYCjwHOBD4nKRtn1EpYnFEHBQRB83X5v1up5kVTMac0scsMpD4tdOOzfnHwqyJmhS/htmiKhf4nghcHJkVZFNH9x1S+8ysTR8SK24qHL/MGqZp8WuYHbINF/hKWkB2ge9lbcvcDRwOkN964IXAyiG20cwKAhoV0AbI8cusYZoWv4Z2AUPFC3w/Dpwv6WayUwQfiog1w2qjmU0XiKnob26zJnL8MmuepsWvoV5RWmH66GrgjcNsU0fJpKqJLzeR4FWpZKypstQ6U+1MJI1Nbi9ZljiiSCV/XZBKDFv+E0wlf02tM5ngtWZi2FTC1VQS18kFqXqlRV22V14vVVZHBGN5JDkKjYlfNmOp5K/JskRC2UEkf32sZqLWx6a2KC+bLL8O+7Gp8rK1iXrJsonypLFPTpaX1dG0+OUpPmaW1KQjTDOzoibFL3fIzKxUPxIrmpmNQtPilztkZlYqS6zYnCNMM7OWpsUvd8jMLEFMNOgI08xso2bFL3fIzKxUBKyfak5AMzNraVr8cofMzEo1bdq4mVlL0+KXO2RmViqAiTG8xYiZWTdNi1/ukA1LIp9YUirvWc2caMlcY6n8ZfNS+cRq5igbQK6xyc1S9VI5wxLrTKTHSeUaq5uHbBBldU01KKCZlZlK5AxLlaVzjZWXPZ0YmambayyVT+zRZD6x8nq/niwvWztRvs5HJxL1EgFz7frydT4x0f8A1qT41ZyWmtnQRYiJmFP6MDMbV/2IX5IWSbpd0gpJp3Uo307S1yXdKGmZpBMLZXdJulnSUklLum3LI2RmltSkazDMzIp6iV+S5gJnAW8AVgHXS7osIm4tLPZHwK0RcYyknYDbJX0pItbl5a+regs1d8jMrFQAE1MeCTOz5ulD/DoYWBERKwEkXQQcCxQ7ZAFsI0nA1sBDwESdjTnSmlmpoLch/27D/fkyh+VD+sskfa/vH8LMZqUK8WuhpCWFx8ltq9gVuKfwelX+XtHngBcBq4GbgT+J2HBxYQBXSbqhw7qfwSNkZlYu6g/5Vxnul7Q9cDawKCLulrRz7402M6NK/FoTEQclyjtVbr9r/JHAUuD1wN7AtyR9PyIeBV4VEavzuPYtSbdFxDVlG/MImZmVag35lz262DDcn19P0RruLzoBuDgi7gaIiPv7/RnMbHbqMX5BdiC5e+H1bmQjYUUnksWwiIgVwJ3AvgARsTr///3AJWQxsZQ7ZGZWqpVYsexBesi/ynD/C4AdJF2dD+u/e7CfyMxmiwrxq5vrgX0kPVfSAuA44LK2Ze4GDgeQtAvwQmClpK0kbZO/vxXwRuCW1MZ8ytLMkibT14qlhvyrDPfPA15OFtC2AK6VdF1E/HzGDTUza9MlfiVFxISkU4ErgbnAeRGxTNIpefk5wMeB8yXdTBbzPhQRayQ9D7gku9afecAFEXFFanvukPWRUglXkxXLe+pSYp2p5K+JZKzJ5K91680rL0smf51fXm8qUZZK4lo7+etm5fUmynM1dtleol4iB2Ldev1ODBs9XENGteH+VWSduseBxyVdAxwAuENmfTWZSOI6Ge3HCRutT9RbX16NpxI3tX5qqvyf3scTiWFTZXWTv/56YsvSskcTiWEfSyR4XZsImGvXlwe3J/ucGLbH+JWvIy4HLm9775zC89Vko1/t9VaSxbLKfMrSzBLE5NSc0kcXVYb7LwVeI2mepC2BQ4Dlff8YZjYL9RS/hs4jZGZWKqB24Koy3B8RyyVdAdwETAHnRkTyOgszsyp6iV+j4A6ZmZWLbNi/dvUuw/356zOAM+pvxcysgx7j17C5Q2ZmpYLeLoo1MxuVpsUvd8jMLKHy9HAzszHTrPjlDpmZJU1NNSegmZkVNSl+uUNmZqUiIBp0hGlm1tK0+OUO2bCkcpTNSfxgUmWpPGR1yxLtjLplqRxlibKpZK6xRD6x+Ym8YIPINZaoN5nKGZaqVzcP2YL+X8E62aAjTJvdpijPGTb1jJzEVcvKrU/8Y7+ORB6yKP8jTpU9NpnI/ZUoG0SusUdTecgSucaeSOUhW9/nRIo0K365Q2ZmpQIx1aBp42ZmLU2LX0NtqaRFkm6XtELSaSXLHCZpqaRlkr43zPaZ2TNF4jGbOH6ZNU+T4tfQRsgkzQXOAt5AdruU6yVdFhG3FpbZHjgbWBQRd0vaeVjtM7MOAqJBQ/6D4vhl1kANi1/DHCE7GFgRESsjYh1wEXBs2zInABdHxN0AEXH/ENtnZh1EqPQxizh+mTVQk+LXMDtkuwL3FF6vyt8regGwg6SrJd0g6d2dViTpZElLJC1ZH08NqLlmFmTTxsses8hA4tcDD04OqLlm1rT4NcyL+jt9+vbTuPOAlwOHA1sA10q6LiJ+Pq1SxGJgMcC2c3Ycx1PBZpuGAMbwSHIEBhK/Djpgc8cvs0FpWPwa5gjZKmD3wuvdgNUdlrkiIh6PiDXANcABQ2qfmXUQU+WPWcTxy6yBeo1f3SbzSNpO0tcl3ZhP5jmxat12w+yQXQ/sI+m5khYAxwGXtS1zKfAaSfMkbQkcAiwfYhvNbBoRU+WPWcTxy6xxeotfhck8RwH7AcdL2q9tsT8Cbo2IA4DDgH+UtKBi3WmGdsoyIiYknQpcCcwFzouIZZJOycvPiYjlkq4AbiLLx3duRNwyrDZWkkqqOmypZLPJxLCJenNTZYl1ppK/JpLbTqUSwya2NzUvUZZK4ppKqppKKJtI1Fo7+WsqEW3N5K9T8/t8Bqxhma4HZZOJX7NYKmnsZCIJwrpIlKWSv06V//OaSv76RCJopMsSCVcTQerxREbstamyASR/fXpdn7skvcevDZN5ACS1JvPcWlgmgG0kCdgaeAiYIDsg61Z3mqEmho2Iy4HL2947p+31GcAZw2yXmSW4QwY4fpk1Ujp+LZS0pPB6cX6NZ0unyTyHtK3jc2Sj5auBbYC3R8SUpCp1p3GmfjNL28QuO5e0Y0Q8OOp2mNkQpOPXmog4KFFeZTLPkcBS4PXA3sC3JH2/Yt1panXI8qG39fnL+yLif9dZj5mNuQB6uFZM0iLgTLLTfOdGxCfbyg8ju/bqzvytiyPir2tvsHt7vgrcL2nbvD3fHdS2zGzEeoxfVJvMcyLwyYgIYIWkO4F9K9adpu4I2bURcSZkR5s112FmDZC4hCapSnb73Pcj4nd7amR1t0XEn+ftOwtwh8xsE1Y3fuU2TOYB7iWbzHNC2zJ3k6W6+b6kXYAXAiuBRyrUnaZuh+xYSVPAle05dsxsE1P/CLPKBbHDtkjSQ8CNZBfemtmmrIcRsiqTeYCPA+dLupnsNOWH8rQ3dKqb2l7dDtm7yPLr/J6kvSPipJrrMbNxFqD6+caqXtT6Skk3kg3nf6Bb0OrRIuAVwKuAnST9S0S8Z4DbM7NR6S1+ZavoMpknIlYDb6xaN6VWhywi7iUbgqu8ITNrIvUyS6nKRa0/BfaMiLWSjga+BuxTt7XdRMQDZHHLsctsk9c1fo2VrolhJb1L0gOSVkl6T/7eoZL+RtINg2+imY3UVOKRz1IqPIpTxrte1BoRj0bE2vz55cB8SQv70WzHLjPrEr/GSpURso8CR5PNgjpV0rfIZhBcCLx/cE2bPZRK1JpM8JooSyRjTSZ/Ta4zkfw1kcQ1UkljEwleI5GoNZkYNlWWWmcqGWviLyWRjzGZULZ28tdEWSwYQJSpf1Fs1wtiJf0G8KuICEkHkx0k9islhWPXLDOZuB/OZOLq7qnEjzz1F7U+ymPb+lTS2ERi2KcSgeGpKA9ETySCzeOJLNRPTCQSvE4k2pIoe3qivJ3r1pfvl4mJ8rLaGpS2p0qHbG1EXA8g6WPAr4AXRMQjg2yYmY2BANW8KLbiBbFvBf5Q0gTwJHBcPn28Hxy7zGazHuLXKFTpkP2GpJOB2/PHKgc0s1mkh+5RhQtiP0eW6XoQHLvMZrtNbITsr4D9gXcALyW7Z9N/Aj8DfhYRFwywfWY2YmpQQGvj2GU2yzUpfnXtkLVdpIuk3ciC3EvJ7mLuoGa2qeo90/XIOHaZzXINi19dO2SSXglc17quIyJWkc2e8rRxs9mgQUeYRY5dZtak+NU17QXwHuCnki6S9N58VpSZzRKaKn+MOccus1muSfGryinLUwAk7Us2zH++pO3I7gF3BfDDiJgcaCvNbHQadIRZ5NhlZk2KX5Uz9UfEbcBtwGckbQG8Dngb8GngoME0r2Hq5gyru84BtCVSZalcY6kcZYmcaKmyqWRus9KidM6wVFkqL1jNsmQ75ydyHyXKkrnG+pyHTA2bNt5Jh9j1Nhy7rCCVo2yy4w0nWmWJPGRRnlNrfSIwpHKUPZ0INk8m8pA9PZnYXqIsmU8skTPs6fXl9VK5xibXVzlpV13T4lfde1l+gezGvEF2tGlmm6oGHWFW8K/AA8B2wAdH3BYzG7QGxa+6HbJrI+JMAEk79rE9ZjZmxvFaix7cHhF/DiDpLHxAabZJa1L8qtshO1bSFHBlRPy8nw0yszESzcrjU8EiSQ8BN5KN8pvZpqph8avSCVtJcyR9pPDWu4A7gN+TdO5AWmZm46FBN+dt1yF2LQKWA68CdpL0L6NpmZkNRY/xS9IiSbdLWiHptA7lH5S0NH/cImlS0rPysrsk3ZyXLem2rUojZBExJekI4O/y1/eS3SzY+XzMNnFNOsJs1yF2PUAWtxy7zGaBXuKXpLnAWcAbyHIYXi/psoi4tbVMRJwBnJEvfwzwvyLiocJqXhcRa6psbyZTGn4m6a8k9XcahJmNt0g8msGxy2y26i1+HQysiIiVEbEOuAg4NrH88cCFdZs6kwC1O3AcsFrSpZI+LultdTdsZg0QzUqsWMKxy2w26j1+7QrcU3i9Kn/vGSRtSXZJxFent4CrJN0g6eRuG5tJHrL/lm90M+DFZPeDOxj4j6rrMLMGas5IWEeOXWazWDp+LWy7tmtx2z1wOyUxK1vjMWTJpounK18VEasl7Qx8S9JtEXFNWWNmPMsyIp4Gfpo/bJSSyV8Tg5+peolkrMl6NZO/xrx6iWhTSWOTCWXLcxISiV2WSvCaWmcqEW2kkr8mykiUaV6fE8PSqJGwJMeuTd9U4l/fdFm59VEeT9Yngsa6RNBIJ40tL3sqkRj26USweWqyvN66qUSC10TS2HWTiQSvk+X7ZTKRGHZqos+JYekav9ZERCo59CqyEfaW3YDVJcseR9vpyohYnf//fkmXkB0IlnbIhnpNRbfZCoXlXpHPVHjrMNtnZm3yaeNlj9nE8cusYXqPX9cD+0h6rqQFZJ2uy9oXym/J9lrg0sJ7W0napvUceCNwS2pjdfOQzViV2QqF5T4FXDmstplZwiYyQtYLxy+zhuohfkXEhKRTyf6e5wLnRcQySafk5efki74ZuCoiHi9U3wW4RNmZpXnABRFxRWp7Q+uQUZitACCpNVvh1rbl/pjsorhXDLFtZlail1OWkhYBZ5IFs3Mj4pMly70CuA54e0R8pf4WB8bxy6yBer3kIiKekSan0BFrvT4fOL/tvZXAATPZ1jBPWXadrSBpV7Ke5rQP207SyZKWSFqyPp7qe0PNLJeaMt5lyL8wqnQUsB9wvKT9SpYb91GlgcSvBx6c7HtDzSzXQ/wahWF2yKrMVvgs8KGISEapiFgcEQdFxEHztXm/2mdmHfQwbbxqDp/WqNL9fW14fw0kfu20Y2JGiJn1rElpe4Z5yrLKbIWDgIvyc64LgaMlTUTE14bSQjN7hi4Xv6amjXcaVTpk2ro3jiq9nvE+zef4ZdZATZp8NMwO2YbZCmS3XToOOKG4QEQ8t/Vc0vnANxzMzEYo6HZRbGra+IxGlZRKqzJ6jl9mTdM9fo2VoXXIZjBbYfQSebOGvr05NfOJ1d1eMrdZKtdYebVUPrFkveQ665Wl8oINoizmJXIfpT5Dop7mlkeYOYl6dYjOvaqKNplRpUbFr03c1AD+hZ1MXFA0mbiyZzIRpKYS9ermIUvlGkuVpXKNpfKJTSTyiU0k6k1NJfJETiYiSqJeHT3Gr6Eb5ghZpdkKhfffO4w2mVlaD9dabFKjSo5fZs0zjteKlRlqh8zMGqhmQPOokpmNnDtkZrZJ6DEjv0eVzGxkGnZHEXfIzCypSUP+ZmZFTYpf7pCZWVqDjjDNzKZpUPxyh8zMykWzjjDNzDZoWPxyh8zM0hp0hGlmNk2D4pc7ZGZWSjTrCNPMrKVp8csdsmEZp2SzKankr4kEr7WTxibamUz+WjdpbOIj1E02W7eMuankr/XK5s7r882qAzTVoENMsxLJ5K81f+Kp5K/rEn/866L8n971iQSvE4kErxOJoDgxlUhSm1hnMvFtIonrZCKhbCTqpcpqaVj8cofMzJKaNG3czKyoSfErMc5gZpYN+Zc9zMzGWa/xS9IiSbdLWiHptA7lH5S0NH/cImlS0rOq1G3nDpmZpUXiYWY2znqIX5LmAmcBRwH7AcdL2m/a6iPOiIgDI+JA4MPA9yLioSp127lDZmblwiNkZtZQvcevg4EVEbEyItYBFwHHJpY/HriwZl13yMysnMiuwSh7mJmNqwrxa6GkJYXHyW2r2BW4p/B6Vf7eM7clbQksAr4607otvqjfzJKaNEvJzKyoS/xaExEHpap3eK9shccAP4yIh2rUBdwhM7OUAPU5k4aZ2VD0Hr9WAbsXXu8GrC5Z9jg2nq6caV3AHbJmG0Cusbr1UvnEkvVSecFSOcqSuc0SZXVzmyXbWX7QkypL5SFLlc1JlaW2V5cHyKwhpuj/hY1TiT/+yY6DIK16iRxeiXWuTwSp1Pbq5iibTOQom0zkBUt9hmQ+sVQw7XceMug1fl0P7CPpucC9ZJ2uE9oXkrQd8FrgnTOtW+QOmZmVa1hiRTOzDXqMXxExIelU4EpgLnBeRCyTdEpefk6+6JuBqyLi8W51U9tzh8zMknzxvpk1Va/xKyIuBy5ve++cttfnA+dXqZviWZZmVqp1L7i608YrJFU8VtJNeVLFJZJePYCPYWazUK/xa9g8QmZm5SKyRw2FxIhvILvA9XpJl0XErYXFvg1cFhEhaX/gy8C+PbbazKyn+DUK7pCZWVIPR5IbEiMCSGolRtzQIYuItYXlt8JTCMysj8ZxJKyMO2RmVi5Ak8k+0kJJSwqvF0fE4vx5p8SIh7SvQNKbgU8AOwO/01uDzcxy3ePXWHGHzMzS0vEslVixUmLEiLgEuETSbwMfB46YaRPNzDpqTn/MHTIzS+th2viMEiNGxDWS9pa0MCLW1N2omVlLk9L2uEPWT3NqTlpNJWodQBLXlKi7veQ6E2W1k7HWSxqbyKuYrlezLJ3AtjxQKDFXW3PKL4qYkyirq4dp410TI0p6PnBHflH/y4AFwIP1W2vWX6lkrJOphKuJwJAqSyVcTSV/Tdcr314k6k2l6iWSuCYTwyZC1CCu92pS2p6hpr2oMAX+HfkU+Jsk/UjSAcNsn5lNpzyxYtkjJSImgFZixOXAl1tJFVuJFYHfA26RtJRsRubbI8ZzWpTjl1mz9BK/RmFoI2QVp8DfCbw2Ih6WdBSwmA4XAZvZEPVw1NotqWJEfAr4VP0tDIfjl1lDeZZlR1WmwP+osPx1ZNecmNmo+NZJLY5fZk3TsPg1zFOWnabA75pY/n3ANzsVSDo5z+q9ZH081ccmmtl05cP9TQp0fTCQ+PXAg5N9bKKZTdes+DXMEbJKU+ABJL2OLKB1vI1KnudoMcC2c3Ycv71qtikZz0u6hm0g8eugAzb3zjUbpAbFr2F2yCpNgc9vn3IucFREeLaV2Sg1LLHiADl+mTVNw+LXME9ZbpgCL2kB2RT4y4oLSNoDuBh4V0T8fIhtM7MykXjMHo5fZk3UoPg1tBGyiJiQ1JoCPxc4rzUFPi8/B/gosCNwtrJ8WBOJLOBmNgSaatA0pQFx/DJrpl7jl6RFwJlkf/fnRsQnOyxzGPBZYD7Z3Utem79/F/AYMEmFeDDUxLAVpsCfBJw0zDbNlOomfx0jdZO/JuulVlk312wqiWvN7Q0iaWw6oWziMCy1OxP15iTK5qa2V0fQqGnjg7QpxK/ZbCoxJDKIn/hkKuFqzYSyU4mgkSxLtGUykcQ1VS+VUDY5+pSs1+fE5D3GryrpbiRtD5wNLIqIuyXt3Laa11W984gz9ZtZKREeITOzRupD/Oqa7obs7iMXR8TdABFxf92NNX+4x8wGJ4DJKH+YmY2r7vFrYSsFTf44uW0NVdLdvADYQdLVkm6Q9O62FlyVv9++7mfwCJmZJalB08bNzIq6xK81Xa7rqpLuZh7wcuBwYAvgWknX5RN7XhURq/PTmN+SdFtEXFO2MXfIzCwhwKcszayReo5fVdLdrCLr2D0OPC7pGuAA4OcRsRqy05iSLiE7BVraIfMpSzMrF2SJFcseZmbjqvf41TXdDXAp8BpJ8yRtSXb/2uWStpK0DYCkrYA3ArekNuYRMjNLalJiRTOzol7iV5V0NxGxXNIVwE1kczrPjYhbJD0PuCRPgTMPuCAirkhtzx0yM0vzSJiZNVWP8atbupv89RnAGW3vrSQ7dVmZO2Q2Xc2T2MkcZXNSCbcSeW7q5hobo7xnyVxjybKa+ctS9eqIgElfQ2bWSSr3V8pk3Zxhydxmfc7hRbovky5L5S8rr9fv8NW0+OUOmZml+aJ+M2uqBsUvX9RvZuUCmIryRxeSFkm6XdIKSad1KH+HpJvyx48kzWiI38ysVI/xa9g8QmZmCQFTk7VqVrntCHAn8NqIeFjSUcBisllKZmY9qh+/RsEdMjMr1zrCrKfrbUci4keF5a8jy/NjZta73uLX0LlDZmZp6WswFkpaUni9OCIW58873XYkNfr1PuCbtdpoZtZJg64hc4fMzBK6JlBM3Xqkym1HsgWl15F1yF49s/aZmZVpVgJrd8jMrFwAk7Wvwahy2xEk7Q+cCxwVEQ/W3ZiZ2TS9xa+hc4fMzBJ6yuOz4bYjwL1ktx05obiApD2Ai4F35TfjNTPrE+chMxtfg0goW9cAkr/2PzEsRNQLaFVuOwJ8FNgRODu/xchE4hSomZVIJYZNJWqtm9y2EXqIX6PgDpmZpfVwhNnttiMRcRJwUu0NmJmleITMzDYJEY2apWRmtkHD4pc7ZGaWFA26KNbMrKhJ8csdMjNLaNa0cTOzjZoVv3wvSzMr15o2XvYwMxtXfYhf3e7Hmy9zmKSlkpZJ+t5M6hZ5hMzMSkVEo4b8zcxaeo1fVe7HK2l74GxgUUTcLWnnqnXbeYTMzJJiKkofZmbjrMf4teF+vBGxDmjdj7foBODiiLgbICLun0HdaRo/QvZYPLTmW+sv/MWo29HBQmDNqBtRgdvZX+Pczj1nWuExHr7yP6e+vDCxyLh+1ka44aan18x99n+NY/yC8f4tF7md/TWu7RxE/No8cS9eqHY/3hcA8yVdDWwDnBkRX6xYd5rGd8giYqdRt6ETSUuakODS7eyvprSzqohYNOo2bMrGNX5Bc37Lbmd/NaWdVfQhflW5H+884OXA4cAWwLWSrqtY9xkrMjMzM7PpqtyPdxWwJiIeBx6XdA1wQMW60/gaMjMzM7Nn2nA/XkkLyO7He1nbMpcCr5E0T9KWZKcll1esO41HyAZncfdFxoLb2V9NaadZN035Lbud/dWUdg5clfvxRsRySVcANwFTwLkRcQtAp7qp7SkalDTNzMzMbFPkU5ZmZmZmI+YOmZmZmdmIuUPWo263RshvqfDr/LYKSyV9dETtPE/S/ZJuKSmXpH/KP8dNkl427Dbm7ejWzpHvT0m7S/qupOX5rTL+pMMyY7E/zVIcv/qrCfErb4dj2DiKCD9qPsgu1LsDeB6wALgR2K9tmcOAb4xBW38beBlwS0n50cA3yXKnHAr8eEzbOfL9CTwbeFn+fBvg5x2+97HYn374UfZw/BpJO8dlfzqGjeHDI2S9mfGtEUYlIq4BHkoscizwxchcB2wv6dnDad1GFdo5chFxX0T8NH/+GNkU513bFhuL/WmW4PjVZ02IX+AYNq7cIetNp1sjtP+oAV4p6UZJ35T04uE0bcaqfpZxMDb7U9JewG8CP24ratL+tNnJ8Ws0xmp/OoaND+ch602VWyP8FNgzItZKOhr4GrDPoBtWw4xv8zAiY7M/JW0NfBV4f0Q82l7coco47k+bvRy/hm+s9qdj2HjxCFlvut4aISIejYi1+fPLyW5CmrrZ6ajM+DYPozAu+1PSfLJA9qWIuLjDIo3YnzarOX4N2TjtT8ew8eMOWW+63hpB0m9IUv78YLJ9/uDQW9rdZcC785k1hwK/joj7Rt2oduOwP/Ptfx5YHhGfLlmsEfvTZjXHryEbl/3pGDaefMqyB1HhtgrAW4E/lDQBPAkcFxFDH/aVdCHZDJ+FklYBfwXML7TzcrJZNSuAJ4ATh93Giu0ch/35KuBdwM2SlubvfQTYo9DOsdifZmUcv0bSzrHYnziGjSXfOsnMzMxsxHzK0szMzGzE3CEzMzMzGzF3yMzMzMxGzB0yMzMzsxFzh8xsxLrdkLjG+iYLNy++rHsNM7N6HL/6x7MszUZM0m8Da8nuG/eSPqxvbURs3XvLzMzSHL/6xyNkZiPW6YbEkvaWdIWkGyR9X9K+I2qemVkpx6/+cYfMNpD0B5J+md/49g5J765QZy9JTxaSC/bahi3yoep1Y3qLlmFZDPxxRLwc+ABw9gzqbi5piaTrJL1pIK0zG0OOYWPD8asGZ+q3ov2B0yPinPy2HpcDX6xQ746IOLAfDYiIJ4EDJd3Vj/U1kbIb/v4W8B/5XVYANsvL3gL8dYdq90bEkfnzPSJitaTnAd+RdHNE3DHodpuNAcewEXP8qs8dMit6KfCV/PmdwLo6K5F0NfAHEXG7pB2B70XESyTtBVwB/AA4FLgR+ALwMWBn4B0R8ZOePsGmYQ7wSKd/IPKbAHe6EXBxmdX5/1fm38VvArMioNms5xg2eo5fNfmUpRW9FLg9v/HsqcCf11zP84H/yp/vD9zcVnZm/v6+wAnAq8mGtT9Sc3ublIh4FLhT0tsguxGwpAOq1JW0g6TW0ehCsnvW3TqwxpqNF8ewEXP8qs8dMgNA0u7ANmRD/PeTBZjza6xnT7Lh56n8rf2BmwqL3BkRN+fly4Bv5zfXvRnYq/YHaLD8hsTXAi+UtErS+4B3AO+TdCPZfjq24upeBCzJ630X+GREzJqAZrOXY9hoOH71j09ZWsv+wDUR8XpJOwC3AK+UdDfwr8BlwKER8fYu6zmQ6cHr5cC/F14/XXg+VXg9xSz9PUbE8SVFi2qs60dkowRms01ZDFtBdprxSuCFwFsLna1ODsQxrDLHr/7xCJm1vBT4GUBEPAxcAPwOcADwtYj4DDBRYT0HAJsDSNqH7Mjo5mQNM7PelcWwVwAXRsSHyUbOduyyHscwGwl3yKxlQzDLfR04miw4XZm/VyWL8IHAnHzI+aPAcuA9/WummVlHZTHsFWQX3wNsFxEPdFnPgTiG2Qg4U78lSToPOAl4FnBaRHygrXwv4ButDM356YHfjIjHetzuXcBBEbGml/WY2eyWX+O0hux04lci4ttt5XvhGGZjYNad77aZiYjfz5+uIZtF1G4S2C5PqvgaYKqXQCZpC7ILROeTXZNhZtaLiYj440S5Y5iNBY+QmZmZmY2YryEzMzMzGzF3yMzMzMxGzB0yMzMzsxFzh8zMzMxsxNwhMzMzMxsxd8jMzMzMRswdMjMzM7MRc4fMzMzMbMT+PxRBEI7syKh5AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 720x216 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# final time\n",
    "plot_concentrations(t[-1])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Input custom particle-size distributions"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to solve the MPM, one must input the area-weighted particle-size distribution $f_{\\text{k},a}$ for each electrode $\\text{k}=\\text{n,p}$ and the minimum and maximum radius limits $R_\\text{k,min}$, $R_\\text{k,max}$. The default distributions $f_{\\text{k},a}$, usable with the Marquis et al. [[6]](#References) parameter set, are lognormals with means equal to the `\"Negative particle radius [m]\"` and `\"Positive particle radius [m]\"` values, and standard deviations equal to 0.3 times the mean.\n",
    "\n",
    "You can input any size distribution $f_{\\text{k},a}(R_\\text{k})$ as a function of $R_\\text{k}$, which we will now demonstrate.\n",
    "\n",
    "Note: $f_{\\text{k},a}(R_\\text{k})$ should ideally integrate to 1 over the specified $R_\\text{k}$ range, although it is automatically normalized within PyBaMM anyway. A distribution such as a lognormal, once restricted to $[R_\\text{k,min},R_\\text{k,max}]$, discretized, and then renormalized, strictly will not integrate to 1 or have the originally desired mean or variance. The mean and variance of the final discretized distribution can be checked as output variables (see below). Having a sufficient number of mesh points in $R_\\text{k}$ or a sufficiently wide interval $[R_\\text{k,min},R_\\text{k,max}]$ should alleviate this issue, however."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameter set (no distribution parameters by default)\n",
    "params = pybamm.ParameterValues(\"Marquis2019\")\n",
    "\n",
    "# Extract the radii values. We will choose these to be the means of our area-weighted distributions\n",
    "R_a_n_dim = params[\"Negative particle radius [m]\"]\n",
    "R_a_p_dim = params[\"Positive particle radius [m]\"]\n",
    "\n",
    "# Standard deviations (dimensional)\n",
    "sd_a_n_dim = 0.2 * R_a_n_dim\n",
    "sd_a_p_dim = 0.6 * R_a_p_dim\n",
    "\n",
    "# Minimum and maximum particle sizes (dimensional)\n",
    "R_min_n = 0\n",
    "R_min_p = 0\n",
    "R_max_n = 2 * R_a_n_dim\n",
    "R_max_p = 3 * R_a_p_dim\n",
    "\n",
    "# Set the area-weighted particle-size distributions.\n",
    "# Choose a lognormal (but any pybamm function could be used)\n",
    "\n",
    "\n",
    "def f_a_dist_n_dim(R):\n",
    "    return pybamm.lognormal(R, R_a_n_dim, sd_a_n_dim)\n",
    "\n",
    "\n",
    "def f_a_dist_p_dim(R):\n",
    "    return pybamm.lognormal(R, R_a_p_dim, sd_a_p_dim)\n",
    "\n",
    "\n",
    "# Note: the only argument must be the particle size R"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# input distribution params to the dictionary\n",
    "distribution_params = {\n",
    "    \"Negative minimum particle radius [m]\": R_min_n,\n",
    "    \"Positive minimum particle radius [m]\": R_min_p,\n",
    "    \"Negative maximum particle radius [m]\": R_max_n,\n",
    "    \"Positive maximum particle radius [m]\": R_max_p,\n",
    "    \"Negative area-weighted \" + \"particle-size distribution [m-1]\": f_a_dist_n_dim,\n",
    "    \"Positive area-weighted \" + \"particle-size distribution [m-1]\": f_a_dist_p_dim,\n",
    "}\n",
    "params.update(distribution_params, check_already_exists=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ef24b8d5c32e4fa19f8665f56d2d1864",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, description='t', max=3511.9522766513314, step=35.11952276651331),…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<pybamm.plotting.quick_plot.QuickPlot at 0x7fd393666a60>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sim = pybamm.Simulation(model, parameter_values=params)\n",
    "sim.solve(t_eval=[0, 3600])\n",
    "\n",
    "sim.plot(output_variables=output_variables)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The discretized size distributions can be plotted as histograms. Only the area-weighted distribution has been input, but the corresponding number and volume-weighted ones are also given as output variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaEAAAEeCAYAAAAjNKpiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA1nUlEQVR4nO3deZwU1bn/8c/DiCIyAioqsnvFhXVQQFxAXFA0KmpcILnGJRH1hhjvjVw1Jopb4s9LEi/GiJgYNIpCUBOSYNwJElEBHQmLXhFRBhCRHQRleX5/nJqxaLpnumd6pmb5vl+veU1X1alTp6q6+ulz6nQdc3dERESS0CjpAoiISMOlICQiIolREBIRkcQoCImISGIUhEREJDEKQiIikhgFoRgzG2tmP024DD82s9/mOc+BZlaSQ/ppZva96PW3zeyFPJZlvpkNjF6PMrPH85h33o9dltt90szOi15fbmYzykn7nJldVmOFS5iZPWNmg8tZ3t/M3q/JMpXHzNqb2SYzKygnzSYzO7SGy7XEzE6ryW1mw8zGm9ldVcljj3wVprYzsyXAQcB2YAewAHgMGOfuOwHc/ZoaLtNA4HF3b1s6z91/VpNlqIi7PwE8UVE6MxsPlLj7TyrIr2s+ylVbjp2Z9QB6At/KJr27n5llvg50dvdFVShebXAP8CDw93QL3f014IgaLVE53P0ToFnptJlNI7zPfhtL0yzNqrVettdoTWtoNaFz3L0Q6EC4OG4EflddGzOzBhPkK1KPj8XVwBNeh371XZPnwt3fAvY1s941tU2pnMSuUXdvEH/AEuC0lHl9gZ1At2h6PHBX9PoA4K/AOmAN8BrQKFrWDngGWAWsBn4dzb8c+Cfwq2idu4C9gNHAJ8BKYCywN7APsCXa/qbo7xBgFOGbF8CvY8s2EWpxo6JlhwBPR2X4CLgutl97R/uyllDjG0n4BpTp2AwC3gPWR9v8B/C92D7NiF5btG+fRWnnAt2A4cA24KuonH+JHfMbo3RfEmreZech2tfJwERgI/A20DNWLgcOi02Pj45phccuSn8uMD86h9OAo1LeDzdEZVsflaFJRec+zbFbDJwYm74cmBGd87XRuTkztnxa7NgeFh3r9cDnwMRo/vRo3zdH+3ZJNP8qYFFUpinAIbF8Twfej/L6TZpzmPq+/DfgFcL793NCbbdFyvEZGR2fzYQvawcBz0Xn6iWgZZS2CfB4lNc6YBZwUCyvh4HbMhy/gcTem+WdlzTrlu7X/VHa94BTY8sPiY7Tmui4XZVy7c8GNhCuy19G8ztGx34P4G5Cq8nW6DyUXucenbt+wKdAQSzf84G50etGwE3Ah9GxmQTsV851eDZQHB3D14Ee6T6/KsoXODFafx2wNDpOuVyj5V03vQjX6cbo3DxF9JlZ0T5k3O8kA0NN/pEmCEXzPwGujV6P5+sg9HNCwGgc/fUnfAgXAO8SLuh9CBfgibGLYjvwg+hk7g3cR7gQ9gMKgb8AP093AUbzRhH7II3NLyIEnF7Rm3AOcCuwJ3Ao4cPwjCjtPYQPzv0IAXNe6nZi+R5AuBAvjPbzP6N9SBeEzoi22yI6FkcBrVOPXcoxL47KsHeai2kU4cIo3fYNhA/txvGLPZZf/PyUe+yAwwkfnoOivP+b8EG0Z6wcbxE+qPYDFgLXlHfu0xy7faIytorNuzzap6sI75VrgeWl67NrEHoSuCU6n2Xvowz7fgohWBxN+GJzPzA95RxeQHjf/TAqQ/wcpr4vD4uOzV5AK0Lguy/l3L1BCDxtCF883ia8//YiBLDborRXE97XTaN9PgbYN5bXfwHPZHj/7XIeyzsvadYt3a//jM7TJYRgtF+0/B+EgNyEr6+fU6NlM4FLo9fNgH7R647Rsd8j9XylOzeEQDAotuyPwE3R6+ujY9g2OmYPAU9m2Jejo2N8bHQML4uOxV5prpuM+QLtCQFiWHRM9geKsr1GKee6if4+jh3vCwnvs7uy2YdMfw2tOS6d5YQ3e6ptQGugg7tvc/fXPBzpvoQLZKS7b3b3re4evxG93N3vd/fthG9QVwH/6e5r3H0j8DNgaC4FNLNWwJ+AH7j7O0AfwgffHe7+lbsvJnzbLM33YuDuaJtLgTHlZH8WsMDdJ7v7NkLQ/DRD2m2EQHok4UN1obuvqKD4Y9x9qbtvybB8TmzbvyR8YPSrIM9sXAL8zd1fjPIeTbjIjk8p23J3X0P4EC2K5mc696laRP83psz/2N0fdvcdwKNRXgelWX8boWn4kDTvo1TfBh5x97fd/UvgZuA4M+tIOIfz3f2Z6H03ht3PYdn70t23uPui6Nh86e6rCMf+pJR17nf3le6+jPCl5k13fyfa/rOEgFS6H/sTPph3uPscd98Qy2dj7FhlI9N5SeczQvDc5u4TCbXBb5hZO0KN4Mbo2BYDvwUujZX5MDM7wN03ufsbOZQv7knCBz5mVkg4F09Gy64GbnH3kuiYjQIuzNDsdRXwkLu/GR3DRwk1k3TXQnn5fht4yd2fjI7J6mjfyxO/Rsu7bvoRgk/p8Z5MqPVWZh/KKAiFb3lr0sz/H8I3gBfMbLGZ3RTNb0f4kNmeIb+lsdetCN8O55jZOjNbR7hB2yrbwplZY0KT1QR3fyqa3QE4pDTPKN8f8/UH3SEp5fi4nE3skjb6sF2aLqG7v0JornsAWGlm48xs3wp2IW1e6ZZ76CBSEpWpqg4htt9R3ksJ57tU/IP6C76+IZ3p3KdaF/0vTJlflq+7fxG9THcz+78JNcq3ol6DV2bcm933ZxOhKaYN6c9ham/IXc6DmR1oZk+Z2TIz20BoTjsgZZ2Vsddb0kyX7tMfgOeBp8xsuZndG71vSxXy9bHKRqbzks6ylC8IHxOOxyFA6Re/+LLS8/9dwrf+98xslpmdnUP54iYAF5jZXoSa6NvuXnqeOgDPxq7RhYTmvXRfSDoAP0q5ptuR/looL992hNpZLuLvjfKum0NIf7wrsw9lGnQQMrM+hIO72zdQd9/o7j9y90OBc4D/MrNTCSekfTk38eIn6HPCxdrV3VtEf83969416b5dp7qf8E0y3qNlKfBRLM8W7l7o7mdFy1cQTn6p9uXkv0taM7OUdXfh7mPc/RigK+EiHlnBvlS0j/FtNyI0MSyPZn1BCOKlDs4h3+WEi6I079L9WlbBeuWd+9R0mwkX/OEV5ZlhO5+6+1Xufgjh2+1vzOywDMlT92cfQu1jGeEcto0ts/h06eZSpn8ezevh7vsC/04IiJXZj23ufru7dyF8Yz4b+E4syVGEJuzq0Cba31LtCcdqObBfVDuJL1sWlfkDdx8GHAj8P2BydExTlfs+c/cFhA/iMwk9JCfEFi8l3A+MX6dNopplqqWE1ot42qbu/mSGtJnyXUq435e2uFnML++6WUH6412ZfSjTIIOQme0bffN5inAP4V9p0pxtZodFB3wD4ZvGDkJ79QrgHjPbx8yamNkJ6bYTfYt4GPiVmR0Y5dvGzM6IkqwE9jez5hnKeTWhieRbUV6l3gI2mNmNZra3mRWYWbcoqEK4UXmzmbU0s7aEewGZ/A3oamYXRIH1Onb9sI+Xp4+ZHRt9y91MaG7cEduXyvx24pjYtq8nVN9Lm0aKgW9F+zeYXZuLyj12hGPwDTM7NSrvj6K8X6+oQOWc+3SmsnszVlbM7KLo/EDoxOBkPp4TgCvMrCj61v0zQvPYEsI57G5m50XH8ftkOIcxhYQb1OvMrA1ff5mozH6cbGbdLfy2ZgOhqSt+vE4idGioDgcC15lZYzO7iBDwpnpohn4d+Hl0jfYg1H6eiMr872bWKrqu1kV5pTvH2byvJxCumwGEe0KlxgJ3m1mHaJutzGxIhjweBq6Jri+LPlu+kRJEs8n3CeA0M7vYzPYws/3NrCiHfSnvuplJuAd3XZT3BYTbE5XZhzINLQj9xcw2EiL2LYR28CsypO1M6AG0iXDwf+Pu06J2/nMIN3Y/ITR7XFLONm8kNO28ETV7vET0uwh3f4/Qfrw4qr6mVluHEd40yy38QG6Tmf04VoYiwo38zwnt3aUfyLcTvp19BLxAaC5Jy90/By4idGZYHe33PzMk35fwRlsb5b+a0GYMofdUl2g//lTO8Uj1Z8LxW0tor78gaouGcIP9HMKHxLcJ98VKy13usXP39wnf7u8nHJ9zCF30v8qiTGnPfYa044Bvp3w7zFYf4E0z20TovPJDd/8oWjYKeDTat4vd/WXgp4QekSsI33aHwi7n8F7COelC6Pn1ZTnbvp1wI3k9IYg9U4nylzqY0GS8gdA09A9C815pa8NmD121q8ObhPP1OaE324XuvjpaNozQ0WA54R7Wbe7+YrRsMDA/Ovb/Cwx1961p8v9fwv2WtWaW6d7qk4QOFq9E5yK+7hRCs+5GwperY9Nl4O6zCfdUfk24FhYROl6kkzFfD79zOosQPNYQvsj1jNar8Bot77qJrp0LonKtJVy3z8TWzWUfypT22BGRSjKzCcAkd/9T0mWBsmbNEuDb7v5qwmV5Gvidu0+thrwvJ/RcOzHfeUvNqa8/IBSpMe6e1dMSqlPUxPsm4R7kSML9ncr2+Mobd/9m0mWQ2q2hNceJ1FfHETpJlDahnOeZu8WL1BpqjhMRkcSoJiQiIonRPaGYAw44wDt27Jh0MURE6ow5c+Z87u5Z/wA/lYJQTMeOHZk9e3bSxRARqTPMrLwnslRIzXEiIpIYBSEREUmMgpCIiCRG94REpNK2bdtGSUkJW7eme+KN1CdNmjShbdu2NG7cuOLEOVAQEpFKKykpobCwkI4dO1K5x+dJXeDurF69mpKSEjp16pTXvNUcJyKVtnXrVvbff38FoHrOzNh///2rpcarICQiVaIA1DBU13lWEBIRkcTonpCI5M1LC1ZWnCgHp3VJNxJ23XLWWWcxYcIEWrRokTHNwIEDGT16NL17995lfnFxMcuXL+ess87KsGZu+dVGCkK1zMZXsh/+pfCUk6uxJCL1x44dOygoKEhk21OnVn4opeLiYmbPnp1zEKpL1BwnInXeeeedxzHHHEPXrl0ZN24cAM2aNePWW2/l2GOPZebMmTz++OP07duXoqIirr76anbsCKN5X3vttfTu3ZuuXbty2223pc3/P/7jP5gyZQoA559/PldeeSUAv/vd7/jJT34CkDH/jh078vnnYcDVO++8kyOPPJJBgwYxbNgwRo8eXbaNP/7xj/Tt25fDDz+c1157ja+++opbb72ViRMnUlRUxMSJE9m8eTNXXnklffr0oVevXvz5z38GYMuWLQwdOpQePXpwySWXsGVL3RnFQ0FIROq8Rx55hDlz5jB79mzGjBnD6tWr2bx5M926dePNN99k//33Z+LEifzzn/+kuLiYgoICnnjiCQDuvvtuZs+ezdy5c/nHP/7B3Llzd8t/wIABvPbaawAsW7aMBQsWADBjxgz69+/PwoULM+Zfavbs2Tz99NO88847PPPMM7s9p3L79u289dZb3Hfffdx+++3sueee3HHHHVxyySUUFxdzySWXcPfdd3PKKacwa9YsXn31VUaOHMnmzZt58MEHadq0KXPnzuWWW25hzpw51XGYq0XiQcjMBpvZ+2a2yMxuSrP8SDObaWZfmtkNsflHmFlx7G+DmV0fLRtlZstiy+pvXVZEGDNmDD179qRfv34sXbqUDz74gIKCAr75zTCw68svv8ycOXPo06cPRUVFvPzyyyxevBiASZMmcfTRR9OrVy/mz59fFmDi+vfvz2uvvcaCBQvo0qULBx10ECtWrGDmzJkcf/zx5eZfasaMGQwZMoS9996bwsJCzjnnnF2WX3DBBQAcc8wxLFmyJO1+vvDCC9xzzz0UFRUxcOBAtm7dyieffML06dP593//dwB69OhBjx49qnQ8a1Ki94TMrAB4ABgElACzzGyKu8ffBWuA64Dz4uu6+/tAUSyfZcCzsSS/cvfRiEi9Nm3aNF566SVmzpxJ06ZNyz6cmzRpUnYfyN257LLL+PnPf77Luh999BGjR49m1qxZtGzZkssvv5ytW7fy5ptvcvXVVwNwxx13cO6557J27Vr+/ve/M2DAANasWcOkSZNo1qwZhYWFGfOPq2gA0b322guAgoICtm/fnjGPp59+miOOOGK3ZXW1q3zSNaG+wCJ3X+zuXwFPAUPiCdz9M3efBWwrJ59TgQ/dvUqPFBeRumf9+vW0bNmSpk2b8t577/HGG2/slubUU09l8uTJfPbZZwCsWbOGjz/+mA0bNrDPPvvQvHlzVq5cyXPPPQfAscceS3FxMcXFxZx77rkAHHfccdx3330MGDCA/v37M3r0aPr3719u/nEnnngif/nLX9i6dSubNm3ib3/7W4X7VlhYyMaNG8umzzjjDO6///6ygPbOO+8AobmwtPlv3rx5aZsUa6uke8e1AZbGpkuAYyuRz1DgyZR5I8zsO8Bs4EfuvjbdimY2HBgO0L59+0psWkRKJdGlevDgwYwdO5YePXpwxBFH0K9fv93SdOnShbvuuovTTz+dnTt30rhxYx544AH69etHr1696Nq1K4ceeignnHBCxu3079+fF154gcMOO4wOHTqwZs2asiCUKf8OHTqUrd+nTx/OPfdcevbsSYcOHejduzfNmzcvd99OPvnksua3m2++mZ/+9Kdcf/319OjRA3enY8eO/PWvf+Xaa6/liiuuoEePHhQVFdG3b99KHs2aZxVVEat142YXAWe4+/ei6UuBvu7+gzRpRwGbUpvYzGxPYDnQ1d1XRvMOAj4HHLgTaO3uV1ZUnt69e3vSg9qpi7bUJQsXLuSoo45Kuhh1xqZNm2jWrBlffPEFAwYMYNy4cRx99NFJFytr6c63mc1x90r/ICnpmlAJ0C423ZYQUHJxJvB2aQACiL82s4eBv1alkCIi+TB8+HAWLFjA1q1bueyyy+pUAKouSQehWUBnM+tE6FgwFPhWjnkMI6Upzsxau/uKaPJ8YF5VCyoiUlUTJkxIugi1TqJByN23m9kI4HmgAHjE3eeb2TXR8rFmdjDhvs6+wM6oG3YXd99gZk0JPeuuTsn6XjMrIjTHLUmzXEREaoGka0K4+1Rgasq8sbHXnxKa6dKt+wWwf5r5l+a5mCIiUg2S7qItIiINmIKQiIgkJvHmOBGpP3L5iUE2kvwZQm0ZDqG+DwWhmpCISJ5leuxOZUydOrXcAFSe4uLiKg0lURMUhESkTluyZAlHHXUUV111FV27duX0009ny5YtDBw4sOxJ1Z9//jkdO3YEYPz48Zx33nmcc845dOrUiV//+tf88pe/pFevXvTr1481a9aU5f34449z/PHH061bN9566y2AjMMpjB8/nosuuohzzjmH008/fZcyaiiIzBSERKTO++CDD/j+97/P/PnzadGiBU8//XS56efNm8eECRN46623uOWWW2jatCnvvPMOxx13HI899lhZus2bN/P666/zm9/8pixwZBpOAWDmzJk8+uijvPLKK7tsT0NBZKZ7QiJS53Xq1ImioiKg/KEQSp188skUFhZSWFhI8+bNy4ZV6N69+y4P/xw2bBgQgsiGDRtYt24dL7zwAlOmTCmrhZQOpwAwaNAg9ttvv922179/f+67776yoSDWrl1bNhTEmDFjePTRR8uGgoBQMznwwAN3ySM+FARQ6aEg0pV9+vTpXHfddUDNDwWhICQidV7pMAgQhkLYsmULe+yxBzt37gTCh22m9I0aNSqbbtSo0S73c1KHRzCzjMMpvPnmm+yzzz5lrzUURHbUHCci9VLHjh3LmpUmT55cqTwmTpwIhFpI8+bNad68ecbhFOI0FET2VBMSkbypTU92v+GGG7j44ov5wx/+wCmnnFKpPFq2bMnxxx/Phg0beOSRRwAyDqdQEQ0FkV6iQznUNhrKQSQ3GsqhZiU9FER9HMpBRESyVB+HglAQqsNUaxJpWOrjUBDqmCAiIolREBIRkcQoCImISGIUhEREJDHqmFAD8v14e5HaatrSaXnNb2C7gXnNb8mSJZx99tnMmzcvr/lWh9mzZ/PYY48xZsyYjGnK25/x48dz+umnc8ghh2S9zSSOj4KQiEgt1Lt37yqN5zN+/Hi6deuWUxBKgprjRKROu/HGG/nNb35TNj1q1Ch+8YtfMHLkSLp160b37t3LHr8TN378eEaMGFE2ffbZZzNt2jQAmjVrxo033sgxxxzDaaedxltvvcXAgQM59NBDy4Zk2LFjByNHjqRPnz706NGDhx56aLdt7Nixg0MPPRR3Z926dTRq1Ijp06cD4QkKixYtyji8wrRp0zj77LMBWLVqFYMGDeLoo4/m6quvpkOHDmXDO+zYsWO3YSwmT57M7Nmz+fa3v01RURFbtmxhzpw5nHTSSRxzzDGcccYZrFixAoA5c+bQs2dPjjvuOB544IGqno6cJR6EzGywmb1vZovM7KY0y480s5lm9qWZ3ZCybImZ/cvMis1sdmz+fmb2opl9EP1vWRP7IiI1b+jQobsEmUmTJnHAAQdQXFzMu+++y0svvcTIkSPLPnSzsXnzZgYOHMicOXMoLCzkJz/5CS+++CLPPvsst956KxDGAmrevDmzZs1i1qxZPPzww3z00Ue75FNQUMDhhx/OggULmDFjBscccwyvvfYaX375JSUlJRx22GHlDg1R6vbbb+eUU07h7bff5vzzzy97ajekH8biwgsvpHfv3jzxxBMUFxezxx578IMf/IDJkyczZ84crrzySm655RYArrjiCsaMGcPMmTNzPvb5kGhznJkVAA8Ag4ASYJaZTXH3BbFka4DrgPMyZHOyu3+eMu8m4GV3vycKbDcBN+a18CJSK/Tq1YvPPvuM5cuXs2rVKlq2bElxcTHDhg2joKCAgw46iJNOOolZs2ZlPUTBnnvuyeDBg4EwvMNee+1F48aN6d69e9kwCS+88AJz584tezjq+vXr+eCDD+jUqdMuefXv35/p06fz0UcfcfPNN/Pwww9z0kknlQ3bUN7QEKVmzJjBs88+C8DgwYNp2fLr79XZDGPx/vvvM2/ePAYNGgSE2lPr1q1Zv34969at46STTgLg0ksv5bnnnsvqGOVL0veE+gKL3H0xgJk9BQwByoKQu38GfGZm38gh3yHAwOj1o8A0FIRE6q0LL7yQyZMn8+mnnzJ06FA+/PDDCteJD/UAuw730Lhx47KhDTIN9eDu3H///Zxxxhm75HvLLbeUPeG6uLiY/v37M3bsWJYvX84dd9zB//zP/zBt2jQGDBhQlk+64RVWrlxZ9rq8Z3ymG8YilbvTtWvX3Wo769atS2wIh1JJN8e1AZbGpkuiedly4AUzm2Nmw2PzD3L3FQDR/wPTrg2Y2XAzm21ms1etWpXDpkWkthg6dChPPfUUkydP5sILL2TAgAFMnDiRHTt2sGrVKqZPn77bk6E7duxIcXExO3fuZOnSpWXDd2frjDPO4MEHH2Tbtm0A/N///R+bN2/m7rvvLhvGAcKwDq+//jqNGjWiSZMmFBUV8dBDD5U9QTuboSFOPPFEJk2aBISa09q1ayssX3wYhyOOOIJVq1aVBaFt27aVNd81b96cGTNmAOw2mmtNSLomlC4E5/JY7xPcfbmZHQi8aGbvufv0XArg7uOAcRCeop3LuiKyq3x3qc5W165d2bhxI23atKF169acf/75zJw5k549e2Jm3HvvvRx88MG7NFWdcMIJdOrUie7du9OtW7ecHwb6ve99jyVLlnD00Ufj7rRq1Yo//elPu6Xba6+9aNeuHf369QNC89yTTz5J9+7dgeyGhrjtttsYNmwYEydO5KSTTqJ169YUFhayadOmjOW7/PLLueaaa9h7772ZOXMmkydP5rrrrmP9+vVs376d66+/nq5du/L73/+eK6+8kqZNm+5Wq6sJiQ7lYGbHAaPc/Yxo+mYAd99teEEzGwVscvfRGfIqW25m7wMD3X2FmbUGprn77kMJpqiuoRxqw++E9ABTqQ4ayqFmfPnllxQUFLDHHnswc+ZMrr322rKaVk2qj0M5zAI6m1knYBkwFPhWNiua2T5AI3ffGL0+HbgjWjwFuAy4J/r/53wXXESkpnzyySdcfPHF7Ny5kz333JOHH3446SLlTaJByN23m9kI4HmgAHjE3eeb2TXR8rFmdjAwG9gX2Glm1wNdgAOAZ6ObansAE9z971HW9wCTzOy7wCfARTW4WyIiedW5c+e094rqg6RrQrj7VGBqyryxsdefAm3TrLoB6Jkhz9XAqXkspohk4O6J97CS6lddt26S7h0nInVYkyZNWL16dbV9QEnt4O6sXr2aJk2a5D3vxGtCIlJ3tW3blpKSEvTzhvqvSZMmtG2brlGqahSERKTSGjduvNsTAkRyoeY4ERFJjIKQiIgkRkFIREQSoyAkIiKJURASEZHEKAiJiEhiFIRERCQxCkIiIpIYBSEREUmMgpCIiCRGQUhERBKjICQiIolREBIRkcQoCImISGIUhEREJDEKQiIikpjEg5CZDTaz981skZndlGb5kWY208y+NLMbYvPbmdmrZrbQzOab2Q9jy0aZ2TIzK47+zqqp/RERkewlOrKqmRUADwCDgBJglplNcfcFsWRrgOuA81JW3w78yN3fNrNCYI6ZvRhb91fuPrp690BERKoi6ZpQX2CRuy9296+Ap4Ah8QTu/pm7zwK2pcxf4e5vR683AguBNjVTbBERyYekg1AbYGlsuoRKBBIz6wj0At6MzR5hZnPN7BEza1mlUoqISLVItDkOsDTzPKcMzJoBTwPXu/uGaPaDwJ1RXncCvwCuzLD+cGA4QPv27XPZdJ2y8ZVXc0pfeMrJ1VQSEZGvJV0TKgHaxabbAsuzXdnMGhMC0BPu/kzpfHdf6e473H0n8DCh2S8tdx/n7r3dvXerVq1y3gEREam8pIPQLKCzmXUysz2BocCUbFY0MwN+Byx091+mLGsdmzwfmJen8oqISB4l2hzn7tvNbATwPFAAPOLu883smmj5WDM7GJgN7AvsNLPrgS5AD+BS4F9mVhxl+WN3nwrca2ZFhOa4JcDVNbZTIiKStaTvCREFjakp88bGXn9KaKZLNYP095Rw90vzWUYREakeSTfHiYhIA6YgJCIiiVEQEhGRxCgIiYhIYhSEREQkMQpCIiKSGAUhERFJjIKQiIgkRkFIREQSoyAkIiKJURASEZHEZPXsODObnmV+W9399CqUR0REGpBsH2DaB7imgjQG/G/ViiMiIg1JtkHodXd/tKJEZvatKpZHREQakKzuCbn7qVmmU1OciIhkrVIdE8zsO/kuiIiINDzlNseZWZd0swkjlT5WLSUSEZEGo6J7Qm8Ak9l9BNMO1VMcERFpSCoKQguBke6+Oj7TzP5WfUUSEZGGoqIgNAjYnDrT3b9RPcUREZGGpNyOCe6+wd13lE6b2YH5LoCZDTaz981skZndlGb5kWY208y+NLMbslnXzPYzsxfN7IPof8t8l1tERKou195xk/O5cTMrAB4AzgS6AMPSdIZYA1wHjM5h3ZuAl929M/ByNC0iIrVMrkEotYNCVfUFFrn7Ynf/CngKGBJP4O6fufssYFsO6w4BSn9c+yhwXp7LLSIieZBrEPI8b78NsDQ2XRLNq+q6B7n7CoDof8ZmRDMbbmazzWz2qlWrsi64iIhUXdJP0U5Xs8o20FVl3a9XcB/n7r3dvXerVq1yXV1ERKog22fHlcp3c1wJ0C423RZYnod1V5pZa3dfYWatgc+qXNIG5qUFK7NOe1qXg6qxJCJSn+VaE7owz9ufBXQ2s05mticwFJiSh3WnAJdFry8D/pzHMouISJ7kVBNy9+y/HmeX33YzGwE8DxQAj7j7fDO7Jlo+1swOBmYD+wI7zex6oIu7b0i3bpT1PcAkM/su8AlwUT7LLSIi+ZFrcxxm1pzQZboX0Cy+rDJP0Xb3qcDUlHljY68/JTS1ZbVuNH81kNWTvxuKd5euy22Fg6ulGCIiu8g5CAF/JNQ8ngW25Lc4IiLSkFQmCPUD9nf31N/tSAOlTgwiUlmV6aI9Azgq3wUREZGGpzI1ocuBqWb2JrDLV2B3vyMfhZLkNZn1etZpt/Y5vhpLIiL1WWWC0N2E3+csIfRYK5XvpymIiEg9V5kgNBQ4vPSxOCIiIpVVmXtCi9n9YaIiIiI5q0xN6A/AFDO7n93vCb2Sl1KJiEiDUJkg9P3o/89S5jtwaNWKIyIiDUnOQcjdO1VHQUREpOHJ+Z5Q9LBQERGRKqtMc9wmM3sPeBcojv4vAW5x9yvyVzSpj/R0BRGJq0wQOhAoiv56AiOA9oRecyIiIlmrzD2hdcC06A8AM7sL0O+GREQkJ/ka3vsu4KY85SUiIg1EZcYTeoBwL6gY+Je7bwVaA1/ktWRSoZzHCBIRqWUqUxNaBpwC/B5YY2YfEHVQMLPzzexIMyvIYxlFRKSeqsw9obIfqZpZY6AL0D36uyr63wpokqcyiohIPZVVTcjM7kw33923ufu77v64u98IzHL3dmhwaBERyUK2NaHrzewRwCpIdx1wW9SDTkREpFzZ3hPaB1iUxV/OTXBmNtjM3jezRWa2Ww87C8ZEy+ea2dHR/CPMrDj2t8HMro+WjTKzZbFlZ+VaLhERqX5Z1YTcPV9duXcRdWB4ABgElACzzGyKuy+IJTsT6Bz9HQs8CBzr7u8TfjBbms8y4NnYer9y99HVUW4REcmPyjwxIZ/6AovcfTGAmT0FDAHiQWgI8Ji7O/CGmbUws9Ypg+qdCnzo7h/XVMGl+ukRPyL1X9JBqA2wNDZdQqjtVJSmDbs+oWEo8GTKeiPM7DvAbOBH7r42XQHMbDgwHKB9+/a5ll9yNG/dzKzTdmtxXDWWRERqg6SDULqODp5Lmuip3ucCN8eWPwjcGaW7E/gFcGW6Arj7OGAcQO/evVO3LQnKJWCdxnnVVxARqTZJB6ESoF1sui2wPMc0ZwJvu3tZ2038tZk9DPw1XwWW3eUSLERE4qqlw0EOZgGdzaxTVKMZCkxJSTMF+E7US64fsD7lftAwUprizKx1bPJ8YF7+iy4iIlWVaE3I3beb2QjgeaAAeMTd55vZNdHyscBU4CxCF/AvgLIxi8ysKaFn3dUpWd9rZkWE5rglaZaLiEgtkHRzHO4+lRBo4vPGxl478P0M634B7J9m/qV5LqbUctOWTss67cB2A6urGCKSo6Sb40REpAFTEBIRkcQoCImISGIUhEREJDGJd0yQ2mku72Wdds8Pss/3q85HVaI0IlJfKQhJvVCc01Dn07JOqZ50ItVLzXEiIpIYBSEREUmMgpCIiCRGQUhERBKjICQiIolR77hK2vjKq0kXQWpALs+kA/WmE8mVakIiIpIYBSEREUmMgpCIiCRG94RqmXdz+uW/iEjdppqQiIgkRjWhBiKXB5LWd7k8Z66oXYtqK4eIKAiJ5JWGGRfJTeLNcWY22MzeN7NFZnZTmuVmZmOi5XPN7OjYsiVm9i8zKzaz2bH5+5nZi2b2QfS/ZU3tj4iIZC/RIGRmBcADwJlAF2CYmXVJSXYm0Dn6Gw48mLL8ZHcvcvfesXk3AS+7e2fg5WhaRERqmaRrQn2BRe6+2N2/Ap4ChqSkGQI85sEbQAsza11BvkOAR6PXjwLn5bHMIiKSJ0kHoTbA0th0STQv2zQOvGBmc8xseCzNQe6+AiD6f2BeSy0iInmRdMcESzPPc0hzgrsvN7MDgRfN7D13n55TAULwGg7Qvn37XFYVEZEqSromVAK0i023BZZnm8bdS/9/BjxLaN4DWFnaZBf9/yxTAdx9nLv3dvferVq1qsKuiIhIrpKuCc0COptZJ2AZMBT4VkqaKcAIM3sKOBZY7+4rzGwfoJG7b4xenw7cEVvnMuCe6P+fq39XRHKj7twiCQchd99uZiOA54EC4BF3n29m10TLxwJTgbOARcAXwBXR6gcBz5oZhP2Y4O5/j5bdA0wys+8CnwAX1dAuSQX2/GBh1mm/6nxUNZZERGqDpGtCuPtUQqCJzxsbe+3A99OstxjomSHP1cCp+S2pNES5PF0B9IQFkVwlfU9IREQaMAUhERFJjIKQiIgkJvF7QiJSMfWkk/pKNSEREUmMgpCIiCRGQUhERBKjICQiIolRx4Q6TEN21z4aOlwkNwpCIvWMetJJXaLmOBERSYyCkIiIJEZBSEREEqMgJCIiiVHHBKm1NPZQ9VMnBkmagpBIQtSdW0TNcSIikiAFIRERSYyCkIiIJEb3hEQkK7l0YgB1ZJDsJF4TMrPBZva+mS0ys5vSLDczGxMtn2tmR0fz25nZq2a20Mzmm9kPY+uMMrNlZlYc/Z1Vk/skIiLZSbQmZGYFwAPAIKAEmGVmU9x9QSzZmUDn6O9Y4MHo/3bgR+7+tpkVAnPM7MXYur9y99E1tS8i1Uk96aS+Srom1BdY5O6L3f0r4ClgSEqaIcBjHrwBtDCz1u6+wt3fBnD3jcBCoE1NFl5ERKom6XtCbYClsekSQi2nojRtgBWlM8ysI9ALeDOWboSZfQeYTagxrU1XADMbDgwHaN++faV2QkR2px/CSjaSDkKWZp7nksbMmgFPA9e7+4Zo9oPAnVG6O4FfAFemK4C7jwPGAfTu3Tt123nxbg5NKSIiDUnSQagEaBebbgsszzaNmTUmBKAn3P2Z0gTuvrL0tZk9DPw1v8WW2kaP+BGpm5K+JzQL6GxmncxsT2AoMCUlzRTgO1EvuX7AendfYWYG/A5Y6O6/jK9gZq1jk+cD86pvF0REpLISrQm5+3YzGwE8DxQAj7j7fDO7Jlo+FpgKnAUsAr4ArohWPwG4FPiXmRVH837s7lOBe82siNActwS4ukZ2SKQWqIs96XT/qOFKujmOKGhMTZk3Nvbage+nWW8G6e8X4e6X5rmYNWYu7yVdBBGRGpN0c5yIiDRgideERERyoaa7+kVBSKQBq4v3j6R+UXOciIgkRjUhaXD0m6KGQ013tZ9qQiIikhjVhEQkK7ncP4K6dw9JtaZkqCYkIiKJUU1IRCRHGmU2fxSERMqRSycGUEeGOHX/lmwoCImIVDPdb8pMQUhEEqdaU8OlICQiUos0tFqTgpBIHumHsNVPtaav1YeApS7aIiKSGNWERKTeUq3pa7W11qQgVAM0UJ2ko6a72kUB62u5/g6qKhSERERyVN8fYVSTFIRE6gDVmuo21bIyUxCqpHdz/CYkUlMUsOq2XGtZ2aqtwS3xIGRmg4H/BQqA37r7PSnLLVp+FvAFcLm7v13euma2HzAR6AgsAS5297U1sT8idYkCVsNRXcGtqhINQmZWADwADAJKgFlmNsXdF8SSnQl0jv6OBR4Ejq1g3ZuAl939HjO7KZq+sab2S6Q+yvU5erlQgGu4kq4J9QUWuftiADN7ChgCxIPQEOAxd3fgDTNrYWatCbWcTOsOAQZG6z8KTCOLILRz40Y2vvJqVgVXjzeR/KnOAFcbKMhmlnQQagMsjU2XEGo7FaVpU8G6B7n7CgB3X2FmB2YqgJkNB4ZHk1/ue+op83LdiTriAODzpAtRjbR/dZv2r+46oiorJx2ELM08zzJNNutWyN3HAeMAzGy2u/fONY+6oD7vG2j/6jrtX91lZrOrsn7Sj+0pAdrFptsCy7NMU966K6MmO6L/n+WxzCIikidJB6FZQGcz62RmewJDgSkpaaYA37GgH7A+amorb90pwGXR68uAP1f3joiISO4SbY5z9+1mNgJ4ntDN+hF3n29m10TLxwJTCd2zFxG6aF9R3rpR1vcAk8zsu8AnwEVZFmlcfvasVqrP+wbav7pO+1d3VWnfLHQ6ExERqXlJN8eJiEgDpiAkIiKJURAiPP7HzN43s0XRExbqFTNbYmb/MrPiqnanrA3M7BEz+8zM5sXm7WdmL5rZB9H/lkmWsSoy7N8oM1sWncNiMzsryTJWlpm1M7NXzWyhmc03sx9G8+vF+Stn/+rL+WtiZm+Z2bvR/t0eza/0+Wvw94Six//8H7HH/wDDUh4dVKeZ2RKgt7vXix/LmdkAYBPhSRrdonn3Amtij2pq6e518lFNGfZvFLDJ3UcnWbaqin4y0drd3zazQmAOcB5wOfXg/JWzfxdTP86fAfu4+yYzawzMAH4IXEAlz59qQrFHB7n7V0Dp43+klnL36cCalNlDCI9oIvp/Xk2WKZ8y7F+94O4rSh9A7O4bgYWEp5/Ui/NXzv7VCx5siiYbR39OFc6fglDmxwLVJw68YGZzoscU1Ue7PKoJyPiopjpshJnNjZrr6mRzVZyZdQR6AW9SD89fyv5BPTl/ZlZgZsWEhwC86O5VOn8KQnl6/E8td4K7H014Ivn3o+YeqVseBP4NKAJWAL9ItDRVZGbNgKeB6919Q9Llybc0+1dvzp+773D3IsJTavqaWbeq5KcglN2jg+o0d18e/f8MeJbQBFnf1OtHNbn7yuji3wk8TB0+h9G9hKeBJ9z9mWh2vTl/6favPp2/Uu6+jjBCwWCqcP4UhLJ7dFCdZWb7RDdIMbN9gNOB+vik8Hr9qKbSCzxyPnX0HEY3tn8HLHT3X8YW1Yvzl2n/6tH5a2VmLaLXewOnAe9RhfPX4HvHAUTdJe/j68f/3J1sifLHzA4l1H4gPKZpQl3fPzN7kjBe1AHASuA24E/AJKA90aOa3L1O3tzPsH8DCU05Thgt+OrSNvi6xMxOBF4D/gXsjGb/mHDfpM6fv3L2bxj14/z1IHQ8KCBUYia5+x1mtj+VPH8KQiIikhg1x4mISGIUhEREJDEKQiIikhgFIRERSYyCkIiIJEZBSEREEqMgJJJG9Jj6gVmkW2Jmp9X0dpNgZm5mm82sWn5nZmavmNlWM5tRHflL7aQgJHVS9OG/xcw2mdlKM/t99Lyuyua1SyBx967uPi0vhc1BdWzXzFpGAWSTmX1hZh+b2XcrmV1Pd78ln+Ur5e6nANdUR95SeykISV12jrs3A44G+gA/yWVlM9ujWkpV+xQBn7t7M3dvCtwMPGRmByRbLBEFIakH3H0Z8BxQOgDcTWb2oZltNLMFZnZ+adqo1nOjmc0FNkePyGkP/CWqKfx3LN1p0et2ZvaMma0ys9Vm9ut05TCzQ8zs6SjdR2Z2XaYyR2VYFpXxfTM7Nc12L4nKVPr3pZlNy3VbhCD0dmz6H4THrlR5OAEzu8XMHoxNtzSzbRZG4FxiZiOj4Qs2m9nvzOwgM3su2u+X6vKQBpIfCkJS55lZO+As4J1o1odAf6A5cDvweMoDJIcB3wBauPswwrOuzolqCvem5F0A/BX4GOhIGGvqqTRlaAT8BXg3SnMqcL2ZnZEm7RHACKCPuxcCZxCeJ7YLd58YlakZcAiwGHgyl21FehFG+CR6+OTPo+lFGdLnojtQHJsuAt53963R9DcJoxYfDpxD+LLwY8Jz8RoB5QVPaQAUhKQu+5OZrSMMMfwP4GcA7v5Hd1/u7jvdfSLwAbs+On+Muy919y1ZbKMvIQCMdPfN7r7V3dPdOO8DtHL3O9z9K3dfTHhk/9A0aXcAewFdzKyxuy9x9w8zFSAKOhOAae7+UI7bghAYfmhmG4C1hAHHBnt+HhyZLgi9G5u+PxrGYBnhwZ5vuvs77v4l4cG6vfJQBqnDGkqbuNRP57n7S6kzzew7wH8Rai4AzQjfvEstTV2nHO2Aj919ewXpOgCHREGxVAHhg3cX7r7IzK4HRgFdzex54L9Kx31K426gkK9rDVlvy8z2Ao4CjnT3D83sm4ShBrZVsD8VioY++TfCE6NL9WTXoLQy9npLmulKdSaR+kM1IalXzKwDoVYwAtjf3VsQxm6Jj6CbWgMor0awFGifRSeGpcBH7t4i9lfo7melS+zuE9z9REJAceD/ZdifoYTmwwvdvTRw5LKtbsCXhKY83P1pQvPjN2PbGGhmfzGzKWY228x6VrCvpboAy9z9iygfIww58W55K4nEKQhJfbMP4UN9FYCZXUHUYaEcK4FDMyx7izAc8z0WBghsYmYnZEi3IepwsLeZFZhZNzPrk5rQzI4ws1OiWspWQo1gR5p0vYD7CTW+VZXZFqG5a15K09tU4NyUdC2BIcClwF0ZjkWq7sCBZvZvFgY4u5MQVJdkub6IgpDUL+6+APgFMJMQXLoD/6xgtZ8DPzGzdWZ2Q0p+Owg31A8j1CBKgEvSbLc0XRHwEfA58FtC54hUewH3RGk+Jdyj+XGadEMIwWFGrIfcczluqwiYmzLv78AgM2sSm/eOBwuBg9Pkk0534HlCZ4NFhOO9GKiW3xFJ/aRB7UQaOAtPaLgTGEDoxTba3c9Jk24roWlvjLv/1MyeA34bNfHloxwvAv2At9z91HzkKbWfOiaICMB6Qrfvg4Gr0iVw9yYps7oDC/NVAHcflK+8pO5QEBIRgPfc/YaKkwXRj0wPJHR/F6k0BSERyZm7rwX2TLocUvfpnpCIiCRGveNERCQxCkIiIpIYBSEREUmMgpCIiCRGQUhERBKjICQiIolREBIRkcQoCImISGL+P7F2vDj5vYiqAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# The discrete sizes or \"bins\" used, and the distributions\n",
    "R_p = sim.solution[\"Positive particle sizes [m]\"].entries[\n",
    "    :, 0\n",
    "]  # const in the current collector direction\n",
    "# The distributions\n",
    "f_a_p = sim.solution[\n",
    "    \"X-averaged positive area-weighted particle-size distribution [m-1]\"\n",
    "].entries[:, 0]\n",
    "f_num_p = sim.solution[\n",
    "    \"X-averaged positive number-based particle-size distribution [m-1]\"\n",
    "].entries[:, 0]\n",
    "f_v_p = sim.solution[\n",
    "    \"X-averaged positive volume-weighted particle-size distribution [m-1]\"\n",
    "].entries[:, 0]\n",
    "\n",
    "\n",
    "# plot\n",
    "width_p = (R_p[-1] - R_p[-2]) / 1e-6\n",
    "plt.bar(\n",
    "    R_p / 1e-6,\n",
    "    f_a_p * 1e-6,\n",
    "    width=width_p,\n",
    "    alpha=0.3,\n",
    "    color=\"tab:blue\",\n",
    "    label=\"area-weighted\",\n",
    ")\n",
    "plt.bar(\n",
    "    R_p / 1e-6,\n",
    "    f_num_p * 1e-6,\n",
    "    width=width_p,\n",
    "    alpha=0.3,\n",
    "    color=\"tab:red\",\n",
    "    label=\"number-weighted\",\n",
    ")\n",
    "plt.bar(\n",
    "    R_p / 1e-6,\n",
    "    f_v_p * 1e-6,\n",
    "    width=width_p,\n",
    "    alpha=0.3,\n",
    "    color=\"tab:green\",\n",
    "    label=\"volume-weighted\",\n",
    ")\n",
    "plt.xlim((0, 30))\n",
    "plt.xlabel(\"Particle size $R_{\\mathrm{p}}$ [$\\mu$m]\", fontsize=12)\n",
    "plt.ylabel(\"[$\\mu$m$^{-1}$]\", fontsize=12)\n",
    "plt.legend(fontsize=10)\n",
    "plt.title(\"Discretized distributions (histograms) in positive electrode\")\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Vary standard deviation as an input parameter\n",
    "You may define the standard deviation (or other distribution parameters except for the min or max radii) of the distribution as a pybamm \"input\" parameter, to quickly change the distribution at the solve stage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define standard deviation in negative electrode to vary\n",
    "sd_a_p_dim = pybamm.Parameter(\n",
    "    \"Positive electrode area-weighted particle-size standard deviation [m]\"\n",
    ")\n",
    "\n",
    "# Set the area-weighted particle-size distribution\n",
    "\n",
    "\n",
    "def f_a_dist_p_dim(R):\n",
    "    return pybamm.lognormal(R, R_a_p_dim, sd_a_p_dim)\n",
    "\n",
    "\n",
    "# input to param dictionary\n",
    "distribution_params = {\n",
    "    \"Positive electrode area-weighted particle-size \"\n",
    "    + \"standard deviation [m]\": \"[input]\",\n",
    "    \"Positive area-weighted \" + \"particle-size distribution [m-1]\": f_a_dist_p_dim,\n",
    "}\n",
    "params.update(distribution_params, check_already_exists=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "78a1747658144f7c86c975c70290a56b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, description='t', max=1.9444444444444444, step=0.01944444444444444…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<pybamm.plotting.quick_plot.QuickPlot at 0x7fd392c3d610>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Experiment with a relaxation period, to see the effect of distribution width\n",
    "experiment = pybamm.Experiment([\"Discharge at 1 C for 3400 s\", \"Rest for 1 hours\"])\n",
    "\n",
    "sim = pybamm.Simulation(model, parameter_values=params, experiment=experiment)\n",
    "solutions = []\n",
    "for sd_a_p in [0.4, 0.6, 0.8]:\n",
    "    solution = sim.solve(\n",
    "        inputs={\n",
    "            \"Positive electrode area-weighted particle-size \"\n",
    "            + \"standard deviation [m]\": sd_a_p * R_a_p_dim\n",
    "        }\n",
    "    )\n",
    "    solutions.append(solution)\n",
    "\n",
    "\n",
    "pybamm.dynamic_plot(\n",
    "    solutions,\n",
    "    output_variables=output_variables,\n",
    "    labels=[\"MPM, sd_a_p=0.4\", \"MPM, sd_a_p=0.6\", \"MPM, sd_a_p=0.8\"],\n",
    ")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Check the distribution statistics\n",
    "The mean and standard deviations of the final discretized distributions can be investigated using the output variables `\"Negative area-weighted mean particle radius\"` and `\"Negative area-weighted particle-size standard deviation\"`, etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The mean of the input lognormal was: 1e-05\n",
      "The means of discretized distributions are:\n",
      "Positive area-weighted mean particle radius [m] 9.972515783613799e-06\n",
      "Positive area-weighted mean particle radius [m] 9.673853099212895e-06\n",
      "Positive area-weighted mean particle radius [m] 9.124186918191047e-06\n"
     ]
    }
   ],
   "source": [
    "print(\"The mean of the input lognormal was:\", R_a_p_dim)\n",
    "print(\"The means of discretized distributions are:\")\n",
    "for solution in solutions:\n",
    "    R = solution[\"Positive area-weighted mean particle radius [m]\"]\n",
    "    print(\"Positive area-weighted mean particle radius [m]\", R.entries[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The standard deviations of the input lognormal were:\n",
      "4.000000000000001e-06\n",
      "6e-06\n",
      "8.000000000000001e-06\n",
      "The standard deviations of discretized distributions are:\n",
      "Positive area-weighted particle-size standard deviation [m] 3.918218937679725e-06\n",
      "Positive area-weighted particle-size standard deviation [m] 5.180362201055076e-06\n",
      "Positive area-weighted particle-size standard deviation [m] 5.815728559306213e-06\n"
     ]
    }
   ],
   "source": [
    "print(\"The standard deviations of the input lognormal were:\")\n",
    "print(0.4 * R_a_p_dim)\n",
    "print(0.6 * R_a_p_dim)\n",
    "print(0.8 * R_a_p_dim)\n",
    "print(\"The standard deviations of discretized distributions are:\")\n",
    "for solution in solutions:\n",
    "    sd = solution[\"Positive area-weighted particle-size standard deviation [m]\"]\n",
    "    print(\"Positive area-weighted particle-size standard deviation [m]\", sd.entries[0])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare to SPM and DFN\n",
    "The MPM can also be easily compared to PyBaMM models with a single particle size. The standard output variables are computed in the MPM, averaging over the particle size domain."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "fbc7fc1fe7434b48bbcb69376b28f58a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, description='t', max=3500.0, step=35.0), Output()), _dom_classes=…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<pybamm.plotting.quick_plot.QuickPlot at 0x7fd3934ec970>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "models = [pybamm.lithium_ion.SPM(), pybamm.lithium_ion.MPM(), pybamm.lithium_ion.DFN()]\n",
    "\n",
    "# solve\n",
    "sims = []\n",
    "for model in models:\n",
    "    sim = pybamm.Simulation(model)\n",
    "    sim.solve(t_eval=[0, 3500])\n",
    "    sims.append(sim)\n",
    "\n",
    "# plot\n",
    "pybamm.dynamic_plot(sims)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model options\n",
    "The MPM is compatible with the current collector and thermal models (except the \"x-full\" thermal option). Currently, the MPM is not compatible with the various degradation submodels in PyBaMM (i.e. SEI models, particle cracking/swelling, or lithium plating)."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fickian diffusion vs Uniform profile\n",
    "One can choose from Fickian diffusion or a uniform concentration profile within the particles. Teh default is \"Fickian diffusion\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1af4a9e7a8c34514829081483c00f571",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, description='t', max=3500.0, step=35.0), Output()), _dom_classes=…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<pybamm.plotting.quick_plot.QuickPlot at 0x7fd3a30ff880>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model_Fickian = pybamm.lithium_ion.MPM(name=\"MPM Fickian\")\n",
    "model_Uniform = pybamm.lithium_ion.MPM(\n",
    "    name=\"MPM Uniform\", options={\"particle\": \"uniform profile\"}\n",
    ")\n",
    "\n",
    "sim_Fickian = pybamm.Simulation(model_Fickian)\n",
    "sim_Uniform = pybamm.Simulation(model_Uniform)\n",
    "\n",
    "sim_Fickian.solve(t_eval=[0, 3500])\n",
    "sim_Uniform.solve(t_eval=[0, 3500])\n",
    "\n",
    "pybamm.dynamic_plot([sim_Fickian, sim_Uniform], output_variables=output_variables)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1D current collector model\n",
    "Add another macroscale dimension \"z\", employing the \"potential pair\" option solving for the potential in the current collectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0f7e454c447545d69e8dae1f8e2b1900",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.0, description='t', max=1.0, step=0.01), Output()), _dom_classes=('w…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<pybamm.plotting.quick_plot.QuickPlot at 0x7fd3307d5190>"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# choose model options\n",
    "model_cc = pybamm.lithium_ion.MPM(\n",
    "    options={\n",
    "        \"current collector\": \"potential pair\",\n",
    "        \"dimensionality\": 1,\n",
    "        \"particle\": \"uniform profile\",  # to reduce computation time\n",
    "    }\n",
    ")\n",
    "\n",
    "# solve\n",
    "sim_cc = pybamm.Simulation(model_cc)\n",
    "sim_cc.solve(t_eval=[0, 3600])\n",
    "\n",
    "# variables to plot\n",
    "output_variables = [\n",
    "    \"X-averaged negative particle surface concentration distribution [mol.m-3]\",\n",
    "    \"X-averaged positive particle surface concentration distribution [mol.m-3]\",\n",
    "    \"X-averaged positive electrode interfacial current density distribution [A.m-2]\",\n",
    "    \"Negative current collector potential [V]\",\n",
    "    \"Positive current collector potential [V]\",\n",
    "    \"Voltage [V]\",\n",
    "]\n",
    "pybamm.dynamic_plot(sim_cc, output_variables=output_variables)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## References\n",
    "\n",
    "The relevant papers for this notebook are:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] Joel A. E. Andersson, Joris Gillis, Greg Horn, James B. Rawlings, and Moritz Diehl. CasADi – A software framework for nonlinear optimization and optimal control. Mathematical Programming Computation, 11(1):1–36, 2019. doi:10.1007/s12532-018-0139-4.\n",
      "[2] Marc Doyle, Thomas F. Fuller, and John Newman. Modeling of galvanostatic charge and discharge of the lithium/polymer/insertion cell. Journal of the Electrochemical society, 140(6):1526–1533, 1993. doi:10.1149/1.2221597.\n",
      "[3] Charles R. Harris, K. Jarrod Millman, Stéfan J. van der Walt, Ralf Gommers, Pauli Virtanen, David Cournapeau, Eric Wieser, Julian Taylor, Sebastian Berg, Nathaniel J. Smith, and others. Array programming with NumPy. Nature, 585(7825):357–362, 2020. doi:10.1038/s41586-020-2649-2.\n",
      "[4] Toby L. Kirk, Jack Evans, Colin P. Please, and S. Jonathan Chapman. Modelling electrode heterogeneity in lithium-ion batteries: unimodal and bimodal particle-size distributions. arXiv:2006.12208, 2020. URL: https://arxiv.org/abs/2006.12208, arXiv:2006.12208.\n",
      "[5] Toby L. Kirk, Colin P. Please, and S. Jon Chapman. Physical modelling of the slow voltage relaxation phenomenon in lithium-ion batteries. Journal of The Electrochemical Society, 168(6):060554, jun 2021. URL: https://doi.org/10.1149/1945-7111/ac0bf7, doi:10.1149/1945-7111/ac0bf7.\n",
      "[6] Scott G. Marquis, Valentin Sulzer, Robert Timms, Colin P. Please, and S. Jon Chapman. An asymptotic derivation of a single particle model with electrolyte. Journal of The Electrochemical Society, 166(15):A3693–A3706, 2019. doi:10.1149/2.0341915jes.\n",
      "[7] Peyman Mohtat, Suhak Lee, Jason B Siegel, and Anna G Stefanopoulou. Towards better estimability of electrode-specific state of health: decoding the cell expansion. Journal of Power Sources, 427:101–111, 2019.\n",
      "[8] Valentin Sulzer, Scott G. Marquis, Robert Timms, Martin Robinson, and S. Jon Chapman. Python Battery Mathematical Modelling (PyBaMM). Journal of Open Research Software, 9(1):14, 2021. doi:10.5334/jors.309.\n",
      "[9] Robert Timms, Scott G Marquis, Valentin Sulzer, Colin P. Please, and S Jonathan Chapman. Asymptotic Reduction of a Lithium-ion Pouch Cell Model. SIAM Journal on Applied Mathematics, 81(3):765–788, 2021. doi:10.1137/20M1336898.\n",
      "[10] Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Evgeni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, and others. SciPy 1.0: fundamental algorithms for scientific computing in Python. Nature Methods, 17(3):261–272, 2020. doi:10.1038/s41592-019-0686-2.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "pybamm.print_citations()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pybamm",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.16"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  },
  "vscode": {
   "interpreter": {
    "hash": "187972e187ab8dfbecfab9e8e194ae6d08262b2d51a54fa40644e3ddb6b5f74c"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
